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ABSTRACT 



M-Layer, the tropospheric propagation effect prediction program by NRaD (formerly 
NOSC), is revised for greater accuracy, speed and stability. This is achieved through 
converting the extended complex number representation into the representation by the 
complex exponent, improving the accuracy in Airy function computation, introducing a 
new mode locating algorithm and implementing a consistency checking procedure for 
determining the proper method to evaluate the height gain function. The revision has 
been documented and the new program source code has been delivered to NRaD. It is 
recommended that the mode search protocol, not just the mode locating algorithm 
introduced in this revision, be completely revised Unlike the current approach of 
blanketing the whole possible region until exhaustion, modes should be searched 
according to their range attenuation rates one by one along a well defined path. This 
should result in a faster and even more stable program. The program size can also be 
reduced. 
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I. INTRODUCTION 



M-Layer is a FORTRAN program for computing the propagation factor of an 
electromagnetic (EM) wave in a stratified atmosphere. It is desirable to extend the 
capability of this program to include a layer of random medium representing the air- 
ocean interface where the thickness of this layer cannot be ignored, where the EM 
propagation and scattering are so strongly coupled that clutter and propagation 
effects within this layer cannot be dealt with separately, and where the grazing angle 
of the EM wave incident into this layer is so small that the curvature of the earth 
cannot be neglected. To achieve this goal, there are many basic theoretical problems 
which have to be answered. First of all, the effect of the earth curvature in this 
program is taken care of through the classical earth-flattening approximation [Ref. 
1], but the result [Ref. 2] does not agree with the more recent diffraction theory of 
Fock [Ref. 3] near the surface of the earth. Then there is the question about the 
better method to model the atmospheric refractive index profile, either piecewise 
linear or quadratic, to be resolved by a new earth-flattening approximation under 
development at NFS. The new approximation will also determine the functions to be 
used for the representation of the EM fields in each layer through uniform 
asymptotic theories. Within some proper region, these new functions are expected to 
reduce to the Airy functions utilized by M-Layer. The evolutionary nature of this 
effort prompted this review to improve the inner workings of the M-Layer program. 
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In particular, the subroutines to search for the modes and those for evaluating the 
Airy functions will remain as an important part of a program investigating questions 
about EM wave propagation by solving the related boundary value problem. 

It can never be overemphasized that a boundary value problem which includes 
a layer of random medium or some range dependent inhomogeneity, set up according 
to the Maxwell equations, will include backscattering in its solution. This is in sharp 
contrast to those numerical procedures based on the parabolic approximation to the 
wave equation for which the backscattering is completely ignored. 

In what follows, the M-Layer program and the reasons for replacing the 
extended complex numbers with their complex exponent representations are 
discussed, together with some other problems encountered and resolved during this 
investigation. 

A. M-LAYER 

In M-Layer, the index of refraction of the atmosphere is assumed to be height 
dependent and is approximated with a continuous piecewise linear profile. The 
classical earth-flattening approximation is utilized to allow the use of the cylindrical 
coordinate system while retaining the effect of the curvature of the earth. This is 
done simply by substituting the index of refraction with the modified index of 
refraction, which also has a piecewise linear profile [Ref. 1]. 

The source of the EM radiation is assumed to be either a vertical electric 
dipole or a vertical "magnetic dipole’, with the latter providing an approximation to 
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the radiation of a horizontal electric dipole. The dipole is located along the positive 
z-axis of the cylindrical coordinate system while the origin is sitting on the ground. 
The x-y plane is the "flattened" earth surface. After carrying out the Hankel 
transform along the radial direction, the resulting spectrum of the Hertzian dipole 
field within each layer of a linear segment of the modified refractive index profile is 
reduced to a linear combination of the Airy functions. Specifically, the layers are 
numbered to increase with height, with the first layer being the one right above the 
ground. The spectrum of the Hertzian dipole field is proportional to the product of 
the values, at the transmitter height and at the receiver height respectively, of the 
height-gain function. At a height within the i-th layer, the height-gain function is 
given by [Ref. 4]; 

where p is the radial component of the propagation vector and is also the spectral 
variable of the Hankel transform; hence it is the same throughout all layers. It is a 
complex variable whose imaginary part represents the radial attenuation rate of the 
spectral component of the Hertzian dipole field. Under the classical earth-flattening 
approximation, the spectrum of the Hertzian dipole field contains a discrete portion 
and a branch cut. The discrete spectrum gives rise to the creeping wave modes 
diffracted by the earth surface and the dielectric waveguide modes supported by the 
layered atmosphere. The contribution from the branch cut is usually negligible, 
especially for the field in the shadow of the earth. The M-Layer program locates the 
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discrete spectrum for modes having a radial attenuation rate below a predetermined 
value. Contributions from these modes determine the propagation factor of the wave. 

The variable q, in the i-th layer is a dimensionless linear function of height z 
with the free space wavenumber k, the modified index of refraction w,- at the lower 
boundary z = z,-, the slope of the modified index of refraction and p as 

parameters: 



<lr 



( v\y p2\ 



k 

a.- , 



(2) 



The height dependence of the field is given in terms of the functions and 

^2(^1')’ which are proportional to the Airy functions Aii—q^e’'’^^^) and Ai{—q^ 

respectively. Of these two functions, at a height so large that (/,• is large and positive, 
represents a downward going wave and €^‘*'^^^k^{q^)+k2{qi) represents an 
upward going wave. The coefficients A- and B- are determined by the conditions on 
the continuity of the Hertzian dipole field and its derivative across layer boundaries 
and by the normalization condition that the integral of the square of the height-gain 
function over all height equals unity. 

To fulfill the radiation condition, the highest layer is given the same refractive 
index as the free space above it and only the outgoing wave is allowed within this 
layer. Below the "flattened" earth surface, the field is assumed to be a plane wave 
propagating downward. Hence only the normalization factors are required in the 

highest layer and in the ground. By assigning B, to unity in the highest layer, all the 
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coefficients and 5 , can be determined, according to the boundary conditions, to 
within a multiplicative factor for This multiplicative factor is then deduced from 
the normalization condition. This procedure can also be carried out from the ground 
level up. That these coefficients can be computed either from the highest level down 
or from the lowest level up is a result of the fact that p belongs to the discrete 
spectrum of the Hertzian dipole field. Consequently, agreement between these two 
ways of evaluating the and Bj coefficients confirms that a mode has been located 
accurately. 

B. EXTENDED COMPLEX NUMBER REPRESENTATION 

The discrete spectrum of the Hertzian dipole field corresponds to the zeroes 
of the modal function which is a determinant whose elements consist of and 
k2(qi) at the layer boundaries. Numerically, the magnitude of this modal function 
causes overflow and underflow problems as k^iq^) or ^2(^1) becomes exponentially 
large or small for complex q- values. In the M-Layer program, to overcome this 
problem, a complex number is written as a scaled number, which is complex, 
multiplied by a scaling factor which is an integer power of e, the base of natural 
logarithm. This integer is chosen so that the greater of the absolute values of the real 
part and the imaginary part of the scaled number lies within e-K A complex number 
written in this form is called an extended complex number. Multiplication of two 
extended complex numbers requires summing the two integer exponents in addition 
to carrying out the regular complex multiplication of the scaled numbers. Addition 
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of two such numbers is achieved through the use of an addition subroutine: the larger 
scaling factor is factored out of both addends before they are combined. The scaling 
factor is adjusted after each addition and after a sequence of multiplications to make 
sure that the resulting scaled number is still within the desired range. Addition is 
troublesome when the two numbers to be added nearly cancel each other. Under this 
circumstance, the scaling factors of the two numbers are identical and both the real 
parts and the imaginary parts of the scaled numbers are almost equal with opposite 
signs. It is clear that the real part and the imaginary part of the sum lose their 
accuracies to different degrees; hence the phase angle may incur substantial error. 
To remedy this situation, interpolation procedures have to be devised. 

As two complex numbers come close to cancel each other, they must be out of 
phase by almost 180 degrees. By factoring out the square root of their product 
instead of the scale factor, the resulting addends become reciprocal to each other, 
both lying within an identical small angle to, and on the same side of, the imaginary 
axis. They are close to the unit circle, but one is on the inside and the other is on the 
outside. Taking out further a phase factor of x/2 after writing the addends in their 
exponential forms, the exponents become small numbers for which Taylor series 
expansion of the exponential function converges rapidly and can be used for 
interpolating the sum to achieve higher accuracy. Note that after the extra phase 
factor of 7 t/ 2 is removed from the addends, it is actually the difference of the 
resulting two reciprocals which is computed. This procedure effectively picks the 
direction on the complex plane along which the addends are almost opposing each 
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other to carry out their cancellation. The resulting sum has a phase angle nearly 
perpendicular to this chosen direction. 

It is evident that the representation of a complex number by its complex 
exponent of base e provides better phase accuracy for addition. A one-to-one 
correspondence can be achieved by restricting the imaginary part of this exponent to 
within -X and x. This will be called the exponential representation or the complex 
exponent representation henceforth. It is convenient for multiplication: adding the 
complex exponents of the two factors will suffice. Conversion of the M-Layer 
program from the extended complex number to the complex exponent representation 
has been carried out. 

C. CONSISTENCY CHECKING 

As better precision is achieved, problems with the mode search procedure and 
the evaluation of the A, and B, coefficients become severe. They are thoroughly 
investigated and resolved. For mode search, although the division of the region of 
interest into "contour rectangles" and further into square "meshes" and the search 
pattern to move around the sides of a "contour rectangle" to find and follow "phase 
lines" into it are kept, the basic assumption of Shellman and Morfitt [Ref. 5] that 
both the real and the imaginary parts of the modal function are linear along every 
edge of a mesh square is completely abandoned. For the evaluation of the A^ and B,- 
coefficients, the "test for evanescence" conditions have been removed. A consistency 
condition to determine whether to evaluate the coefficients from the ground level up 
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or from the top level down has been fomulated and incorporated into the program. 
This accomplishment leads to the relaxation of mode locating accuracy requirement 
which, combined with the improved precision of the revised program, makes the first 
order Newton-Raphson iteration unnecessary. The specific changes in the program 
and the resulting gains in speed, accuracy and execution stability are discussed in the 
following chapters. Recommendation to completely revise the mode search protocol 
to do without the "contour rectangles" is also provided. 
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II. PROGRAM REVISIONS 



M-Layer is structured into three parts: setup, mode search and propagation 
factor evaluation. The main input is the modified refractive index values at specified 
heights so that a piecewise linear profile can be constructed. If the mode locations 
for the particular profile are available from a previous run of the program, they can 
also be included in the input and the mode search procedures will be bypassed. The 
various ranges and transmitter and receiver heights for which propagation factors are 
desired are also specified. The subroutine WVGSTDIN is called to input the 
information from an ASCII data file. The program then computes the constants to 
be used for mode search and propagation factor evaluation. The mode search is 
performed with the subroutine FNDMOD. The MODSUM subroutine is then 
invoked to first compute the /!,• and S,- coefficients as explained in the Introduction, 
then compute the propagation factor and the propagation loss. The complete 
program structure is given in Figures 1 and 2. There are several other subroutines 
which are not included in these and other figures, such as DHORIZ for computing 
the horizon distance between a transmitter and a receiver for reference purpose; 
CHKMOD, a maintenance routine for removing zero from reported mode locations 
by older versions of the program; or A02H20, a routine to compute the atmospheric 
absorption coefficient due to oxygen and water vapor. They will not be discussed as 
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Figure 1 Original M-layer subroutines structure 






Figure 2 Original M-layer subroutines structure (continued) 
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they do not contribute directly to the main purpose of this program of locating the 
modes and computing the propagation factor. 

The program structure has been altered as shown in Figures 3 and 4. Since the 
and B- coefficients have to be evaluated only once, they are now obtained through 
a call to the subroutine ABCOEF directly from the main program right after the 
modes are located. Several subroutines are dropped in this revision for various 
reasons; The subroutines NORME and NORMRE are eliminated because they are 
no longer needed due to the change in complex number representation; The 
subroutines NOMSHX, FDFDTX and DXDETR are not used because the modes 
are now located with adequate precision without further iteration; The subroutine 
ADDX is not listed separately because it is called only once and has been reduced 
to only a few lines which are placed where the subroutine is called in the original 
program. On the other hand, changes in the mode search algorithm require the 
addition of two new subroutines: SURFO is a modified and simpler version of SURF; 
ROOTS replaces QUAD. Due to the change in complex number representation, all 
subroutines listed below FNDMOD and MODSUM have been revised, including 
their input/output lists. But except for SURFO and ROOTS, the utilities of these 
subroutines are the same as those of the original ones. Descriptions of these 
subroutines can be found in the report by Yeoh [Ref. 4]. 

The most significant changes have been made in XCADD, XCDAIT and 
XCDAIG for adopting the complex exponent representation and improving 
computation speed and accuracy; in FZEROX, FINDFX, ROOTS and SURFO for 
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Figure 3 New M-layer subroutines structure 
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Figure 4 New M-layer subroutines structure (continued) 
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stabilizing and simplifying the mode search algorithm; and in ABCOEF for 
implementing the criteria to determine the reliable manner for evaluating the /!,• and 
coefficients. These changes are discussed in the sections below. The source code 
listings of the completely new subroutines XCADD and ROOTS and the significantly 
revised subroutines FZEROX and ABCOEF, which are compiled with Microsoft 
FORTRAN version 5. 00, are attached as Appendices A through D. Validation of the 
revised program has been carried out at 9.6 GHz for all the 21 profiles listed in 
Yeoh [Ref. 4]. 

A. ADDITION SUBROUTINE 

XCADD is the subroutine implementing the addition of complex numbers 
under the representation by their exponents. Given the double precision complex 
numbers Zj and as the exponents of the addends, this subroutine returns the 
exponent of the sum. Since a double precision number has an accuracy of 53 bits, if 
the real parts of Zj and Z 2 differ by more than 53 bits, the exponent of their sum will 
simply be the one of the greater real part. When cancellation becomes serious, the 
square root of the addends is factored out first. Then the four-term Taylor series 
expansions of the resulting reciprocals are summed up. Since the leading term of the 
sum of the Taylor series is a good estimate of the sum of the reciprocals and the 
relative error of the four-term Taylor series sum is proportional to the fourth order 
of this leading term, the threshold for invoking this interpolation procedure is set at 
the highest possible value of 2“*^^ allowed under double precision. Experimenting 
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with this procedure shows that this interpolation improves accuracy as long as the 
threshold is set at a number between 2“^'^ and 

B. AIRY FUNCTION EVALUATION 

Similar to the original program, the evaluation of the Airy function adopted the 
algorithm prescribed by Schulten, et. al. [Ref. 6], In the new program, changes are 
made to follow the advice of Schulten, et. al. concerning the region within which 
Taylor series expansion, instead of the faster Gaussian quadrature, has to be used to 
achieve double precision accuracy. Other changes in implementing the algorithm are 
described below. 

1. XCDAIT 

Due to the similarity in their Taylor series coefficients, the Airy function 
and its derivative are evaluated within a single loop. The relative accuracy of the 
derivative of the Airy function is set at the double precision limit of 

2 . XCDAIG 

Six term Gaussion quadrature is used for evaluating the Airy function and 
its derivative outside the circle of radius 4.97 centered at (0.90, 2.80) on the complex 
plane. The use of four-term quadrature outside a radius of 15 from the origin 
suggested by Schulten, et. al. is not adopted. The six-term quadrature in this range 
retains a higher accuracy while overall speed improvement by using both the four- 
term and the six-term quadrature appears to be minimal. 
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C. MODE LOCATING 



As explained in the Introduction, the modes are located at the zeroes of the 
modal function. These zeroes are located on the upper complex plane. Here qjj 
is the value of q] on the earth surface, which, according to Eq.(2) of Chapter 1, is a 
linear function of p^. For a horizontally propagating mode, p/k is close to unity. The 
maximum range attenuation rate specified for the desired modes, which corresponds 
to a limit on the imaginary part of p, determines approximately the upper bound for 
the imaginary part of the qjj complex plane to be searched for modes. The Shellman 
and Morffit mode search procedure first divides the search region horizontally into 
"contour rectangles" each of which spans 160 meshes along the real q^ direction. A 
mesh is a square whose size is an adjustable parameter of the order 10“"* at 9.6 GHz 
for most of the cases considered herein. This parameter is determined by the 
frequency and the slope of the modified index of reflection in the lowest layer of the 
profile. The search commences at the top left comer of the "contour rectangle" whose 
left edge has a real coordinate value close to the difference of the real parts of the 
qii values with the minimum modified index of refraction and the index near the 
surface substituted into Eq.(2) of Chapter 1. After the search over the initial 
rectangle is completed, the program moves to search the next rectangle until a 
specified maximum number of modes are found or a specified number of "contour 
rectangles" have been searched. 

The search for zeroes makes use of the fact that a real function changes sign 
when it crosses a simple zero. Since a zero of a complex valued function F(q) is 
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where both its real part and imaginary part vanish, a necessary condition for a point 
q^j to be a zero is that it is on the intersection of two curves defined by Im{F(q)}=0 
and Re{F(q)}=0. The program searches around a "contour rectangle" for a sign 
change in Im{F(q)} across an edge of a mesh bordering the side of the "contour 
rectangle" to determine that a line of Im{F(q)}=0 has been encountered. The search 
then follows this line into the meshes within the "contour rectangle', checking each 
mesh to see if a curve Re{F(q)}=0 enters the mesh under investigation. All these 
steps make use only of the assumption that the zeroes of the modal function are 
simple. Once both the curve Im(F(q)}=0 and the curve Re{F(q)}=0 are determined 
to be present within a mesh, the location of their possible interception is estimated. 
An algorithm for this estimate is required. 

Shellman and Morffit [Ref. 5] introduced a further assumption that the 
functions Re{F(q)} and Im{F(q)} are both linear along the edges of a mesh. Based 
on this assumption, they try to estimate the locations where the curve Im{F(q)}=0 
enters and leaves a mesh square and the location of qj^j if a curve Re{F(q)} =0 also 
enters the same mesh. It is obvious that information about the locations where the 
curves enter and leave the mesh square is not essential. Furthermore, in the 18 m 
duct height case, the scheme causes the search path to loop around four contiguous 
meshes until the search is broken up by the limit on the number of meshes to be 
investigated. Replacing their technique requires major changes in the subroutines 
involved. A new subroutine ROOTS is provided to estimate the location of the 
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intersection of the curves Im{F(q)}=0 and Re{F(q)}=0. These changes eliminate 
the looping problem. 

Another problem is encountered in the 40 m duct height case when a large 
number of zeroes are found in the lower half complex qjj plane. These zeroes 
appear to belong to the reflection coefficient on the wrong sheet of the branch cut 
and are not waveguide modes. This happens because the search region has been 
extended below the real q^ axis to avoid the singularity in SURF. The problem with 
this singularity should have been solved within SURF, especially because it occurs 
only when the derivative of the subroutine output variable gamma with respect to q^ 
is computed. Since this derivative is not needed during mode search, the extension 
of the search region to the negative qjj plane is unnecessary. A simplified routine, 
SURFO, is introduced which is exactly the same as SURF except that it does not 
evaluate the derivative of gamma. By using this subroutine instead of SURF, the 
search path in the revised program does not avoid the real and the imaginary axes. 

1. FNDMOD 

The search region is limited to the upper half qjj plane. All the modes 
found are ordered according to their range attenuation rates before those numbered 
beyond the maximum modes allowed are abandoned. 

2. FZEROX 

Since the curve Im(F(q))=0 enters into a mesh square through an edge, 
the values of Im{F(q)} must change sign over the end points of either one or all 
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three other edges. When there is only one other edge across which Im{F(q)} changes 
sign at its end points, it is the edge across which the curve Im{F(q)}=0 exits the 
mesh square. Ambiguity arises when all edges indicate a change of sign at their end 
points. When this occurs, a "right turn rule" is adopted which assumes that the curve 
exits the edge to the right of the one along which it enters the mesh square. Such a 
rule avoids the retracing of the search path when the mesh square is revisited as 
entering this same mesh square from the left side of an edge after exiting from its 
right side requires a crossing of the Im{F(q)}=0 curve, which is prohibited under the 
simple zero assumption. On the other hand, the actual curve may have turned left 
and then returns to this mesh square, i.e., following a "left turn rule." Under such a 
scenario, this wrong choice would have left a segment of the curve not searched. This 
difficulty has not been observed during testing. In fact the ambiguous situation 
seldom occurs. Note also that, as remarked above, two lines of Im{F(q)}=0 do not 
cross each other unless a higher order zero is present. Hence only a right turn rule 
or a left turn rule for the curve to exit the mesh is allowed. Exiting the opposite edge 
demands a pair of crossing Im{F(q)}=0 curves within the mesh square. This violates 
the assumption that all zeroes are simple. Also note that, the possibility of vanishing 
Re{F(q)} or Im{F(q)} values at the comers of a mesh square is eliminated through 
a small adjustment in FINDFX. 

3. FINDFX 

Both the vertical shift away from the real q^ axis and the horizontal 
offset away from the imaginary axis are unnecessary and have been removed from 
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this routine. Furthermore, as a result of converting to the complex exponent 
representation, the sine and cosine of the argument of the modal function are 
examined for sign changes in FZEROX. This is implemented in FINDFX by 
including the cosine and sine values of the argument of the modal function in the 
output list. To avoid the indeterminate case when either the real or the imaginary 
part of the modal function becomes zero at any comer of a mesh square, the 
argument for computing the cosine and sine values is increased by 2“^^ when this 
occurs. This is equivalent to a consistent small distortion of the particular comer of 
the mesh square. This will not cause any error in locating the zero because FINDFX 
still returns separately the unmodified exponent of the value of the modal function. 

4. ROOTS 

Assuming that the modal function is analytic within the mesh, this 
subroutine utilizes the values of the modal function at the four comers of the mesh 
square to determine the Taylor series expansion coefficients of the modal function 
to the third order. The roots of this cubic polynomial are then located using Cardan's 
solution by radicals. If the higher order coefficients fall below machine resolution for 
a root within the mesh square, these coefficients are regarded as zero and the order 
of the polynomial is reduced and can be solved more expediently. If the function is 
determined to be constant over the mesh square, the center of the square is taken 
as the root location. 
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D. EVALUATINGA, AND 



As discussed in the Introduction, the A, and B,- coefficients can be evaluated 
either from the top level down or from the lowest level up. These two procedures are 
simply called "integration down" and "integration up" respectively in the original 
documentation [Ref. 4]. The location of a mode has been called an eigenvalue. That 
the results of integration down and integration up agree is a manifestation that the 
eigenvalue is located accurately. 

The subroutine ABCOEF evaluates the coefficients A, and B,- for each mode. 
If the range attenuation rate for a mode is greater than 0.1 dB/km, the coefficients 
are evaluated from the lowest layer up. Otherwise, it is evaluated from the top layer 
down. It is obvious that such a rule must be implemented because the results of 
integration up and integration down do not agree for many modes. Efforts are made 
to determine the cause of this discrepancy and to devise a means to resolve it. 

Investigation reveals that inadequate precision in the location of the modes is 
one source of the problem. Since the B, coefficients depend on the A,- coefficients 
while the A,- coefficients are obtained directly, only the A,- coefficients need to be 
examined. The Ay coefficients of the six modes of lowest range attenuation rates for 
all 21 profiles except the one without evaporation duct are computed using 
eigenvalues of different accuracy controlled by the first order Newton-Raphson 
iteration method. Table 1 shows the A, coefficient computed with the new program. 
They are arranged from the top layer down. In the i-th layer, the A,- coefficient 
computed by integration downward depends only on A,.^y in the layer above while 
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TABLE 1. IMPROVING Aj ACCURACYWITH EIGENVALUE (18 M DUCT) 



mode 


4 q-eigenvalue: 


.1888574325176803D*»-00 


.10806787448105980- 


■01 






eigenvalue difference: 


.OOD+OO 


.000+00 


-.490-13 - 


.660-15 


.120-11 - 


.120-10 


.150-06 


.60D-07 


layer U 


A i /down 


A i /down 


A i /down 


A i /down 


A i /down 


layer # 


Ai/up 




Ai/up 


Ai/up 


Ai/up 


Ai/up 


18 


.0261 


.6719 


.0261 


.6719 


.0261 


.6719 


.0261 


.6719 


.0261 


.6719 


18 


.0261 


.6719 


.0261 


.6719 


.0261 


.6719 


.0261 


.6719 


.0261 


.6719 


17 


-.0625 


.6368 


-.0625 


.6368 


-.0625 


.6368 


-.0625 


.6368 


-.0625 


.6368 


17 


-.0625 


.6368 


-.0625 


.6368 


-.0625 


.6368 


-.0625 


.6368 


-.0625 


.6368 


16 


.0139 


.7440 


.0139 


.7440 


.0139 


.7440 


.0139 


.7440 


.0139 


.7440 


16 


.0139 


.7440 


.0139 


.7440 


.0139 


.7440 


.0139 


.7440 


.0139 


.7440 


15 


.1216 


.6353 


.1216 


.6353 


.1216 


.6353 


.1216 


.6353 


.1216 


.6353 


15 


.1216 


.6353 


.1216 


.6353 


.1216 


.6353 


.1216 


.6353 


.1216 


.6353 


14 


.0166 


.5471 


.0166 


.5471 


.0166 


.5471 


.0166 


.5471 


.0166 


.5471 


14 


.0166 


.5471 


.0166 


.5471 


.0166 


.5471 


.0166 


.5471 


.0166 


.5471 


13 


-.1565 


.5310 


-.1565 


.5310 


-.1565 


.5310 


-.1565 


.5310 


-.1565 


.5310 


13 


-.1565 


.5310 


-.1565 


.5310 


-.1565 


.5310 


-.1565 


.5310 


-.1565 


.5310 


12 


-.3842 


.5659 


-.3842 


.5659 


-.3842 


.5659 


-.3842 


.5659 


-.3843 


.5659 


12 


-.3842 


.5659 


-.3842 


.5659 


-.3842 


.5659 


-.3842 


.5659 


-.3842 


.5659 


11 


2.2002 - 


.8081 


-2.2002 


-.8081 


-2.2002 


-.8081 


-2.2002 


-.8081 


-2.1909 


-.8068 


11 


2.2002 - 


.8081 


-2.2002 


-.8081 


-2.2002 


-.8081 


-2.2002 


-.8081 


-2.2002 


-.8081 


10 


5.4648 


.2423 


-5.4648 


.2423 


-5.4648 


.2423 


-5.4654 


.2423 


-4.1810 


-.2161 


10 


5.4648 


.2423 


-5.4648 


.2423 


-5.4648 


.2423 


-5.4648 


.2423 


-5.4647 


.2423 


9 


3.6974 - 


.6979 


-3.6974 


-.6979 


-3.6974 


-.6980 


-3.6783 


-.7012 


-6.4611 


-.2121 


9 


3.6978 - 


.6978 


-3.6978 


-.6978 


-3.6978 


-.6978 


-3.6978 


-.6978 


-3.6977 


-.6978 


8 


.3459 - 


.7982 


.3459 


-.7982 


.3460 


-.7982 


.3482 


-.7926 


-1.9078 


-.9148 


8 


.3459 - 


.7983 


.3459 


-.7983 


.3459 


-.7983 


.3459 


-.7983 


.3459 


-.7983 


7 


.4098 


.8794 


.4098 


.8794 


.4098 


.8794 


.4136 


.8836 


-1.0899 


.5364 


7 


.4097 


.8793 


.4097 


.8793 


.4097 


.8793 


.4097 


.8793 


.4097 


.8793 


6 


.3480 


.8161 


.3480 


.8161 


.3480 


.8161 


.3526 


.8205 


-.5879 


.4005 


6 


.3479 


.8160 


.3479 


.8160 


.3479 


.8160 


.3479 


.8160 


.3479 


.8160 


5 


.2923 


.8304 


.2923 


.8304 


.2923 


.8304 


.2972 


.8358 


-.3490 


.3749 


5 


.2922 


.8303 


.2922 


.8303 


.2922 


.8303 


.2922 


.8303 


.2922 


.8303 


4 


.2359 


.8619 


.2359 


.8619 


.2360 


.8619 


.2408 


.8690 


-.2058 


.3731 


4 


,2358 


.8618 


.2358 


.8618 


.2358 


.8618 


.2358 


.8618 


.2358 


.8618 


3 


.1831 


.8910 


.1831 


.8910 


.1832 


.8910 


.1878 


.9003 


-.1250 


.3753 


3 


.1831 


.8908 


.1831 


.8908 


.1831 


.8908 


.1831 


.8908 


.1831 


.8908 


2 


.1300 


.9149 


.1300 


.9149 


.1301 


.9149 


.1342 


.9275 


-.0734 


.3750 


2 


.1300 


.9146 


.1300 


.9146 


.1300 


.9146 


.1300 


.9146 


.1300 


.9146 


1 


.0586 


.9335 


.0586 


.9335 


.0588 


.9335 


.0618 


.9545 


-.0318 


.3670 


1 


.0586 


.9331 


.0586 


.9331 


.0586 


.9331 


.0586 


.9331 


.0586 


.9331 
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that computed by integration upward depends only on A^_j in the layer below. Hence 
in each layer, the coefficient obtained by integration downward is listed above that 
obtained by integration upward. There are five sets of values listed, with the 
magnitudes given in powers of 10 and the phase given as a multiple of t. They are 
obtained from eigenvalues of decreasing accuracy, the one used to compute the left 
most column being the most accurate. The first set is computed using an eigenvalue 
having a relative accuracy The second set uses an eigenvalue with a relative 

accuracy of 2”^^; The relative accuracy of the eigenvalue for the third set is 2“^^; 
For the fourth set, the first order Newton-Raphson iteration of the mode location is 
set at an absolute accuracy of 0.03 of the mesh size, same as that specified in the 
original program; The eigenvalue for the right most set is the mode location 
estimated by ROOTS without modification by the Newton-Raphson iteration. It is 
clear that, for this mode, the difference between these two methods of computing the 
coefficients becomes negligible as the accuracy in mode location increases. For 
example, in the 8-th layer, the magnitude of A^ computed by integrating downward 
changes from —1.9078 to 0.3482 to 0.3460 to 0.3459, which agrees with the result 
computed by integrating upward. The phase follows the same trend to an agreement 
within O.OOlx.Table 2 shows a similar set of output, but the coefficients fail to agree 
even when the relative accuracy is increased to 2“'^®. Note that the actual difference 
in both the real part and the imaginary part of the two most accurate eigenvalues is 
about 2““*^. Double precision accuracy appears to be insufficient for the coefficients 
computed with these two methods to agree for all modes. Some interesting features 
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TABLE2. IMPROVING A, ACCURACYWITH EIGENVALUE (36 M DUCT) 



mode 


3 q- eigenvalue 


: .31480001647813920+00 


.1479622940572007D 


•02 






eigenvalue difference: 


.380-14 


.360-14 


.380-14 


.360-14 


-.160-09 - 


.220-09 


-.530-07 


.280-07 


layer # A i /down 


A i /down 


A i /down 


A i /down 


A i /down 


layer # Ai/up 




Ai/up 


Ai/up 


Ai/up 


Ai/up 


27 


-.0009 


.6663 


-.0009 


.6663 


-.0009 


.6663 


-.0009 


.6663 


-.0009 


.6663 


27 


.2353 


.7582 


.2353 


.7582 


.2353 


.7582 


.2353 


.7582 


.2353 


.7582 


26 


.0007 


.6678 


.0007 


.6678 


.0007 


.6678 


.0007 


.6678 


.0007 


.6678 


26 


-.0111 


.3659 


-.0111 


.3659 


-.0111 


.3659 


-.0111 


.3659 


-.0111 


.3659 


25 


.0022 


.6657 


.0022 


.6657 


.0022 


.6657 


.0022 


.6657 


.0022 


.6657 


25 


-1.8851 


.3913 


-1.8851 


.3913 


-1.8851 


.3913 


-1.8851 


.3913 


-1.8852 


.3913 


24 


.0001 


.6809 


.0001 


.6809 


.0001 


.6809 


.0001 


.6809 


.0001 


.6809 


24 


-7.4914 


.6081 


-7.4914 


.6081 


-7.4914 


.6081 


-7.4914 


.6081 


•7.4914 


.6081 


23 


-2.9495 


.5951 


-2.9495 


.5951 


-2.9495 


.5951 


-2.9495 


.5951 


-2.9495 


.5951 


23 


-14.5340 


.7973 


-14.5340 


.7973 


-14.5340 


.7973 


-14.5340 


.7973 


-14.5340 


.7973 


22 


-12.1956 


.9278 


-12.1956 


.9278 


-12.1956 


.9278 


-12.1956 


.9278 


-12.1956 


.9278 


22 


-23.5827 - 


.9407 


-23.5827 


-.9406 


-23.5827 


-.9406 


-23.5827 


-.9406 


-23.5827 


-.9406 


21 


-35.2395 - 


.2502 


-35.2395 


-.2502 


-35.2395 


-.2502 


-35.2395 


-.2502 


-35.2396 


-.2501 


21 


•44.4517 - 


.8199 


-45.8691 


.8599 


-45.8691 


.8599 


•47.4590 


-.1252 


-47.4594 


-.1251 


20 - 


131.3304 - 


.9570 


-131.3304 


-.9570 


-131.3304 


-.9570 


-131.3304 


-.9570 


-131.3307 


-.9569 


20 • 


129.0146 - 


.2961 


-127.6070 


.0248 


-127.6070 


.0248 


-122.9124 


-.9081 


-120.9305 


.8279 


19 


-25.6088 - 


.9230 


-25.6088 


-.9230 


-25.6088 


-.9230 


-25.6088 


-.9230 


-25.6088 


-.9230 


19 


-25.6090 - 


.9228 


-25.6184 


-.9241 


-25.6184 


-.9241 


-22.5644 


-.8054 


-20.2166 


.7391 


18 


-13.6970 


.6510 


-13.6970 


.6510 


-13.6970 


.6510 


-13-6970 


.6510 


-13.6970 


.6510 


18 


-13.6970 


.6510 


-13.6970 


.6510 


-13.6970 


.6510 


-13.0618 


.7675 


-10.8148 


.3440 


17 


-7.0384 


.4145 


-7.0384 


.4145 


-7.0384 


.4145 


-7.0384 


.4145 


-7.0384 


.4145 


17 


-7.0384 


.4145 


-7.0384 


.4145 


-7.0384 


.4145 


-7.0308 


.4179 


-6.3129 


.1800 


16 


-3.3146 


.2991 


-3.3146 


.2991 


-3.3146 


.2991 


-3.3146 


.2991 


-3.3146 


.2991 


16 


•3.3146 


.2991 


-3.3146 


.2991 


-3.3146 


.2991 


•3.3146 


.2991 


•3.3116 


.2970 


15 


-2.3132 


.2632 


-2.3132 


.2632 


-2.3132 


.2632 


-2.3132 


.2632 


-2.3132 


.2632 


15 


•2.3132 


.2632 


•2.3132 


.2632 


-2.3132 


.2632 


-2.3132 


.2632 


•2.3127 


.2629 


14 


-1.5669 


.2415 


-1.5669 


.2415 


-1.5669 


.2415 


•1.5669 


.2415 


-1.5669 


.2415 


14 


-1.5669 


.2415 


-1.5669 


.2415 


-1.5669 


.2415 


-1.5669 


.2415 


-1.5668 


.2415 


13 


•1.0838 


.2352 


-1.0838 


.2352 


-1.0838 


.2352 


-1.0838 


.2352 


•1.0838 


.2352 


13 


-1.0838 


.2352 


-1.0838 


.2352 


-1.0838 


.2352 


-1.0838 


.2352 


-1.0838 


.2352 


12 


-.6983 


.2432 


-.6983 


.2432 


-.6983 


.2432 


-.6983 


.2432 


-.6983 


.2432 


12 


-.6983 


.2432 


-.6983 


.2432 


-.6983 


.2432 


-.6983 


.2432 


-.6983 


.2432 


11 


-.3754 


.2712 


-.3754 


.2712 


-.3754 


.2712 


-.3754 


.2712 


-.3754 


.2712 


11 


-.3754 


.2712 


-.3754 


.2712 


-.3754 


.2712 


-.3754 


.2712 


-.3754 


.2712 


10 


-.0102 


.3619 


-.0102 


.3619 


-.0102 


.3619 


-.0102 


.3619 


-.0102 


.3619 


10 


-.0102 


.3619 


-.0102 


.3619 


-.0102 


.3619 


-.0102 


.3619 


-.0102 


.3619 
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can be observed in both tables, which are present in all 120 sets of values computed. 
When disagreement is present in one set of coefficients such as those in either 
Table 1 or Table 2, the change toward smaller differences with improving eigenvalue 
accuracy occurs mainly in one way of computation, but not both. For example, in 
Table l,the values of integration downward improve with better eigenvalue accuracy, 
while those computed by integrating upward change little. In Table 2, the results of 
integration downward are the ones that are holding steady as the accuracy in 
eigenvalue improves. Furthermore, when disagreement occurs, the layer in which the 
coefficient has the smallest magnitude, i.e.,the one having the most negative 
power of 10, divides the table into two parts. The results of two different ways of 
computation agree in the layers above this one if they disagree in those below it, and 
vise versa. No explanation will be attempted. Instead, practical rules are drawn up 
to take advantage of these facts. In Table 1 , the process of integration upward goes 
through the troublesome 10-th layer and produces results which agree with the results 
of downward integration before the downward process goes through the 10-th layer. 
On the other hand, the downward integration is tripped up going across the 10-th 
layer and produces results which fail to agree with the results from upward 
integration. It is clear that the results from upward integration are the correct ones. 
This conclusion is further supported by the fact that improving the accuracy of the 
eigenvalue does not change significantly the results of upward integration. Similar 
argument leads to the conclusion that in Table 2, the results of downward integration 
are the correct values. 
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It can be concluded from the above observations that one of the methods of 



computing the /!,• coefficients converges to the correct value much faster then the 
other. It is also found that this method of faster convergence is always able to arrive 
at the correct values for A- for all the cases under investigation. 

Table 3 lists the statistics of the method of integration which yields the correct 
Ai coefficients for each of the 120 modes investigated. The differences in magnitudes 
and phases in the lowest layer and in the layer below the highest are also listed. 
Since for most of the cases when disagreement in A- values occurs, the correct 
integration is upward, this is used as the default. To decide that downward 
integration should be utilized, the following steps are taken: The first /!,• value of 
downward integration is computed and compared to the value from upward 
integration. If the magnitudes in dB disagree by less than 0.02 dB, their phases will 
be checked. If the phases differ by less than the agreement is deemed 

acceptable and the /I, and coefficients computed from the lowest layer up are 
used. Otherwise, the coefficients are re-evaluated again from the highest layer down. 

Once the correct method of evaluating the A^ and B- coefficients is used, the 
accuracy of the mode location becomes less critical. For all the cases investigated, 
the A( coefficients obtained from mode locations estimated with or without the 
Newton -Raph son first order iteration differ only by 0.06 dB in magnitude and 
0.0013xin phase at most. In fact, few cases show differences more than 0.002 dB and 
0.000 lx. The Newton-Raphson iteration is not needed. Hence the subroutines 
NOMSHX, FDFDTX and DXDETR are removed. 
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TABLES. STATISTICS FOR EVALUATING A, COEFnCIENT 



Duct 


Mode # 


Evaluating Method 


AKI 


m 


Aarg(A^)ln 


height 


























Layer 


Layer 






up 


down 


bottom 


top-1 


bottom 


top-1 




1 


X 














2 


X 














3 


X 












02 


1 


X 




0.172 




0.093 






5 


X 














K 


X 




8362 




13234 






1 


X 














2 


X 














3 


X 




0.008 




0.0002 




04 


1 


X 




1.030 




1.8717 






5 


X 




7.814 




13948 






6 


X 




0.002 




0.0001 






1 


X 














2 


X 




0.002 




0.0004 






3 


X 




0.522 




0.0158 




06 


4 


X 














5 


X 




13J78 




0.4377 






6 


X 




0.002 




0.0001 






1 


X 














2 


X 




0.002 










3 


X 




0.002 




0.0001 


0.0001 


08 


4 


X 




0.016 




0.0026 






5 


X 




4.066 




0.6355 






.. 6 


X 




3.978 




0.6186 





26 



TABLE 3. CONTINUED 1 



Duct 


Mode # 


Evaluating Method 


A 1.4,1 


m 


Aarg(/l.)/7t 


height 


























Layer 


Layer 






up 


dov\Ti 


bottom 


top-1 


bottom 


top-1 




1 


X 














2 


X 








0.0002 






3 


X 








•.0002 




10 


4 


X 




0.04 




0.0008 






5 


4 




0206 




0.0402 


0.0001 




6 


X 




0.002 




0.0001 






1 


X 














2 


X 




0.006 




0.0003 






3 


X 




0.004 








12 


4 


X 




1.808 




0.5661 






6 


X 




1.732 




0.5429 






6 


X 




1.472 




0.0414 






1 


X 














2 


X 




0.002 




0.0001 






3 


X 




0.178 




0.0052 




14 


4 


X 




0.024 




0.0005 






5 


X 




0.004 




0.0001 






6 


X 




0.85 




0.4711 






1 


X 














2 


X 




0.006 




0.0002 






3 


X 




0.004 








16 


4 


X 




0.006 ' 




0.0001 






5 


X 




0.002 




0.0001 






... 6 


X 




0.004 




0.0077 
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TABLE 3. CONTINUED 2 



Duct 


Mode # 


Evaluating Method i 


A 1-1,1 


m 


Aarg(A^)lTz 


height 




































Layer 


Layer 






up 


do\%Tl 


bottom 


top-1 


bottom 


top-1 




1 




X 




0.008 




0.0001 




2 


6 




0.002 




0.0001 






3 


X 








0.0001 




18 


4 


X 














5 


X 




0.016 




0.0003 






6 


6 




0.002 










1 




X 




0.078 




0.0164 




2 


4 














3 


X 




0.002 




0.0001 




20 


4 


6 








0.0008 






5 


X 




0.16 




0.0195 






6 


4 




0.002 




0.0011 






4 




X 




8.«08 




0339 




2 


6 














4 


X 




0.004 








22 


3 


6 




0.016 










5 


X 




0.002 




0.0001 






6 


X 




031 




0.0117 






1 


X 














2 




X 




0.868 




03842 




3 


X 




0.006 




0.0009 




24 


4 


X 




0.002 




0.0001 






5 


X 




0.026 




0.0009 






6 


X 




0.008 




0.0001 
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TABLES. CONTINUED 3 



Duct 


Mode # 


Evaluating Method 


A 1^,1 


m 


Aarg(/l,)/7t 


height 




































Layer 


Layer 






up 


do>^Tl 


bottom 


top-1 


bottom 


top-1 




1 


X 




0.002 


0.002 


0.0001 


0.0001 




2 




X 




4308 




0.121 




3 


X 




0.006 








26 


1 


X 




0.002 




0.0001 






5 


X 








0.0001 






Z 


X 




0.034 




0.0039 






1 




X 




0.028 




0.0014 




2 




X 




4.806 




0.0728 




3 


X 












28 


4 


X 














5 


X 




0.008 


0.002 


0.0002 






1 


X 




0.004 




0.0019 






1 




X 




1.562 




0.0165 




2 


X 














3 




X 




0.718 




03455 


30 


4 


X 














5 


X 




0.004 










$ 


X 




0.724 




0.0522 






1 




X 




3.194 




0.1648 




2 


X 




0.002 










3 




X 




13.12 




0.1026 


32 


4 


X 




0.002 










5 


X 




0J82 




0.0099 






6 


X 




0.002 




0.0001 
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TABLES. CONTINUED 4 











A \A,\ 


m 


Aarg(A^)/n 


Duct 


Mode # 


Evaluating Method 










height 
























Layer 


Layer 






up 


dov^Ti 


bottom 


top-1 


bottom 


top-1 

1 




1 


X 




0.002 


0.002 








2 




X 




13.456 




0.0311 




3 




X 




1.014 




0.2347 


34 


4 


X 














5 


X 




0.03 




0.0006 






6 


X 




0.014 




0.0006 






1 




X 






0.0001 


0.0014 




2 




X 




1.686 




0.2224 




3 




X 




4.724 




( 

0.0919 


36 


4 


X 










\ 




5 


X 




0.006 




0.0001 






6 


X 




0.02 




0.0001 






1 




X 




0.996 




0.0115 




2 




X 




4.974 




0.0152 




3 


X 












38 


4 




X 




5.052 




0.0417 




5 


X 








0.0001 


- i 




6 


X 




0.002 






T 




1 


X 




0.002 


0.002 




f 




2 




X 




3.85 




0.1226 




3 




X 




3.568 




0.1555 , 


40 


4 




X 




3.448 




0.1678 




5 


X 








0.0001 






6 


X 
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III. CONCLUSION AND RECOMMENDATION 



A. Performance 

This revision of M-Layer converts the extended complex number representation 
of an exponentially large or small number into the direct representation by its 
complex exponent. The accuracy of the computation has been improved in two ways; 
First, an interpolation algorithm has been devised when severe cancellation of the 
addends is detected. Secondly, accuracy for the evaluation of the Airy function has 
been improved, not just by summing the Taylor series to double precision resolution 
and by adopting six-term Gaussian quadrature, but also by expanding the region 
within which the more expedient Gaussian quadrature is excluded in favor of the 
more accurate but time-consuming Taylor series summation. The improvement in 
accuracy is most easily seen from Table 1. 

As discussed in the Introduction, evaluating the Aj and 5, coefficients either 
from the lowest layer up (integration up) or from the top layer down (integration 
down) must result in the same values. This property provides a consistency check for 
the accuracy of the computation. For the six modes of lowest range attenuation rates 
of the 20 profiles of different duct heights. Table 1 lists the maximum difference for 
each mode which shows a discrepancy between these two methods of evaluating the 
/!,• coefficients. For each profile, the maximum value in magnitude difference in dB 
among all the layers is listed if it is greater than 2. If the phases of the coefficients 
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TABLE 1. MAXIMUM DIFFERENCE IN A, COEFFICIENT BETWEEN 
INTEGRATION UP AND DOWN 



Duct 

height 

(m) 


Mode 


Difference in A,- coefficient 


Magnitude difference in 
(dB) 


Phase difference over 

0.1 T 


original 


revised 


original 


revised 


02 


3 


5.22 




Yes 




6 


61.16 




Yes 




04 


2 


22.46 


2.3 






5 


106.9 




Yes 




06 


3 


8.62 




Yes 




5 


32.36 








04 


5 


77.84 




Yes 




6 


44.9 




Yes 




10 


5 






Yes 




12 


3 


69.38 




Yes 




5 


46.32 




Yes 




6 


7.46 




Yes 




14 


6 


30.6 




Yes 




22 


1 


8.64 




Yes 




14 


2 


80.48 




Yes 




14 


2 


110.68 




Yes 




28 


2 


lt6.9 


67.68 


Yes 


Yes 


30 


3 


173.28 


143.42 


Yes 


Yes 


32 


1 


11.38 




Yes 




3 


525.04 


188.04 


Yes 


Yes 
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TABLE 1. CONTINUED 



Duct 

height 

(m) 


Mode 

ft 


Difference in coefficient 


Magnitude difference in 
(dB) 


Phase difference over 

0.1 V 


original 


revised 


original 


revised 


33 


2 


37.98 




Yes 




3 


715.7 


209.94 


Yes 


Yes 


36 


2 


112.74 




Yes 




3 


957.92 


231.68 


Yes 


Yes 


38 


2 


107.44 


52.26 


Yes 


Yes 


4 


1249 


255.8 


Yes 


Yes 


40 


3 


167 


112.72 


Yes 


Yes 


4 


823.56 


258.18 


Yes 


Yes 




Magnitude difference within 2dB are not listed. 



deviate more than 0. Ixin any layer, that particular mode is also singled out. The 
location of the mode of the revised program is within a relative accuracy of 
achieved through first order Newton-Raphson iteration. Even though discrepancies 
still exist when the duct is 28 meters or higher, it is clear that the revised program 
computes more accurately than the original one. 

For the cases where the two methods of evaluating the and Bj coefficients 
disagree, it has been observed that one of the methods always leads to Aj values 
which are little changed when the accuracy in mode location is varied, while the 
other method produces A^ values which shift toward the results of the other method 
as the accuracy of mode location improves. Based on this observation, a consistency 
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check is implemented into the program to identify the method which converges 
better. For the 120 cases investigated, when this method of faster convergence is 
used, the A^ coefficients obtained from mode locations estimated with or without the 
Newton-Raphson first order iteration differ only by 0.06 dB in magnitude and 
O.OOlSirin phase at most. In fact, few cases show differences more than 0.002dB and 
0.0001 7T. This allowed the Newton-Raphson iteration to be removed in this revision. 

Table 2 compares the performance between the original and the revised 
programs. The time spent to find the modes has been reduced by an average of 
22.58%. The revised program can always produce the modes found by the original 
program. Moreover, the mode search is stable for the new program: the time it 
requires to search for the modes is about the same for similar profiles. The sudden 
jumps in mode search time for the 24 m and the 40 m cases, which indicate troubles 
during the search, no longer happen. 

With the proper method of evaluating the A- and B, coefficients determined by 
the consistency check, the output of the revised program differs from the original 
program in some cases. The most serious deviation has been observed for the 38 m 
duct height case as shown in Tables 3 and 4. For example, at a range of 36.5 km with 
the transmitter at a height of 25 m and the receiver at 10 m, the coherent path loss 
is 175.93 dB from the original program, and is 167.90dB from the revised program. 
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TABLE 2. OVERALLMODE SEARCH PERFORMANCE COMPARISON 



DUCT 

HEIGHT 

(meters) 


ORIGINAL PROGRAM 


REVISED PROGRAM 


Time 

Improvement 


Time 


Modes 


Time 


Modes 


(X) 


0:00:37 


3 


0:00:35 


3 


5.40% 


02 


0:32:14 


9 


0:31:55 


9 


0.98% 


04 


1:14:12 


25 


1:05:04 


25 


12.31% 


08 


2:10:18 


53 


1:56:50 


53 


10.33% 


08 


0:35:58 


39 


0:29:25 


39 


18.21% 


16 


0:53:24 


59 


0:48:32 


61 


9.11% 


12 


1:09:40 


86 


1:01:44 


89 


11.39% 


14 


1:20:42 


94 


1:11:13 


97 


11.75% 


16 


1:54:35 


95 


1:18:07 


97 


31.82% 


16 


1:45:09 


103 


1:27:15 


104 


17.02% 


20 


1:46:19 


103 


1:34:20 


105 


11.27% 


2t 


1:52:54 


lOS 


1:35:18 


108 


15.59% 


24 


3:42:59 


106 


1:46:47 


107 


52.11% 


26 


2:07:42 


106 


1:43:55 


108 


18.62% 


28 


2:00:05 


107 


1:44:59 


109 


12.57% 


30 


1:59:59 


107 


1:46:19 


108 


11.39% 


•2 


1:55:29 


108 


1:42:58 


110 


10.84% 


34 


2:29:57 


109 


2:15:58 


111 


9.32% 


36 


2:31:40 


109 


2:17:20 


112 


9.45% 


1 30 


2:38:44 


110 


2:18:09 


111 


12.97% 


1 36 


5:41:17 


95 


2:39:39 


111 


53.22% 


Total 


40:23:54 




31:16:22 




22.58% 
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TABLES. ORIGINAL PROGRAM OUTPUT: 38 M DUCT 



frequency = 


9600.0000 mhz 










range 


zt 


zr 


coherent 


incoherent 


coherent 


incoherent 


horizon 


(1cm) 


(m) 


(m) 


mode sum 


mode sum 


path loss 


path loss 


(km) 








(db) 


(db) 


(db) 


(db) 




27.3 


25.0 


4.0 


-15.30 


-15.62 


156.10 


156.43 


28.9 


27.3 


25.0 


6.0 


.62 


-2.35 


140.18 


143.16 


30.7 


27.3 


25.0 


8.0 


-1.11 


-4.21 


141.92 


145.01 


32.3 


27.3 


25.0 


10.0 


-27.26 


-12.66 


168.06 


153.46 


33.6 


36.5 


25.0 


4.0 


-16.94 


-16.62 


160.28 


159.96 


28.9 


36.5 


25.0 


6.0 


-.73 


-2.05 


144.07 


145.39 


30.7 


36.5 


25.0 


8.0 


-2.21 


-3.72 


145.55 


147.06 


32.3 


36.5 


25.0 


10.0 


-32.59 


-14.29 


175.93 


157.64 


33.6 


45.8 


25.0 


4.0 


-19.89 


-16.96 


165.20 


162.26 


28.9 


45.8 


25.0 


6.0 


-2.81 


-1.89 


148.11 


147.19 


30.7 


45.8 


25.0 


8.0 


-4.11 


-3.43 


149.41 


148.74 


32.3 


45.8 


25.0 


10.0 


-28.57 


-15.22 


173.88 


160.52 


33.6 



TABLE 4. REVISED PROGRAM OUTPUT; 38 M DUCT 



frequency = 


9600. OOOC 


1 mhz 










range 


zt 


zr 


coherent 


incoherent 


coherent 


incoherent 


horizon 


(km) 


(m) 


(m) 


mode sum 


mode sum 


path loss 


path loss 


(km) 








(db) 


(db) 


(db) 


(db) 




27.3 


25.0 


4.0 


-14.38 


-15.66 


155.18 


156.47 


28.9 


27.3 


25.0 


6.0 


.42 


-2.37 


140.39 


143.18 


30.7 


27.3 


25.0 


8.0 


-1.52 


-4.21 


142.33 


145.02 


32.3 


27.3 


25.0 


10.0 


-21.20 


-12.51 


162.01 


153.31 


33.6 


36.5 


25.0 


4.0 


-17.32 


-16.60 


160.66 


159.94 


28.9 


36.5 


25.0 


6.0 


-.48 


-2.08 


143.82 


145.42 


30.7 


36.5 


25.0 


8.0 


-1.62 


-3.73 


144.96 


147.07 


32.3 


36.5 


25.0 


10.0 


-24.56 


-14.04 


167.90 


157.38 


33.6 


45.8 


25.0 


4.0 


-20.26 


-16.93 


165.57 


162.23 


28.9 


45.8 


25.0 


6.0 


-3.14 


-1.93 


148.44 


147.23 


30.7 


45.8 


25.0 


8.0 


-4.62 


-3.46 


149.92 


148.76 


32.3 


45.8 


25.0 


10.0 


-25.40 


-14.90 


170.71 


160.21 


33.6 
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B. Recommendation 



The mode search protocol of this program needs to be revised. Since the search 
is limited by the maximum range attenuation rate accepted, it is logical to begin with 
locating the mode of the lowest or the highest attenuation, then proceed to look for 
the next one in the order of increasing or decreasing attenuation rate. Furthermore, 
under the assumption of analyticity over the search region, there should be only one 
connected "phase line" of vanishing real part of the modal function on which all the 
modes are located. The partition of the search region into rectangles as has been 
done in this program tends to cut the "phase line" into segments before the program 
starts to search for the end points of these segments and then follow the segments 
in different directions. It is clear that a better way is to search for one end of the 
"phase line" along a line of a constant attenuation rate in the search region, either 
at the maximum accepted or the minimum possible attenuation, then follow this 
"phase line" all the way to the other end. This technique works even if the "phase 
line" branches off into several directions at a Stokes' point. 
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APPENDIX A: SUBROUTINE XCADD 



This Appendix lists the addition subroutine XCADD which returns the complex 
exponent of the sum when the complex exponents of the addends are given. This is 
a complete re-write of the original subroutine of the same name. 
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1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 

21 

22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 



subroutine xcaddCzx, zlx, z2x) 
c 

c Given zlx and z2x, this subroutine adds the two complex numbers 
c z1=exp(z1x) and z2=exp(z2x) for z=exp(zx) and returns zx. 

c 

c inputs... 

c z1x=complex exponent of the complex number z1 

c z2x=complex exponent of the complex ntinber z2 

c 

c outputs... 

c zx=complex exponent of the complex number z 

c 

c subroutines called... 

c 

implicit real*8 (a‘h,o-z) 

complex* 16 zx, zlx, z2x, zt lx, zt2x, cl ogzh,dsum,czero, cerrx, cone, chpi 
parameter (pi=3. 141592653589793238462643d0, twopi =2. d0*pi , 

+ hpi=0.5d0*pi ,zero=0.d0,c16=1 .d0/6.d0, 

♦ bit 14=1. d0/16384.d0,bit24=bit 14/1024. d0,ctol=bit 14, 

♦ dpi =2259. d0/4294967296.d0/4294 967296. d0,hdpi=dpi/2.d0, 

+ e2m54= -3 . 742994 775023 7048 19d1 , e2p27= * 0 . 5d0*e2m54 , 

+ chpi = (0.d0,1.5 70796326794 8966 1 923 1 32d0 ) , c one= ( 1 . dO , 0 . dO ) , 

+ czero=(0.d0,0.d0),cerrx=( -3. 742994 775023704819d1,0.d0)) 

c cerrx=e2m54=-54*log(2)=exponent below machine accuracy 

dimension ztmp(2) , stmp(2) 
equivalence (ztmp, clogzh) , (stmp,dsum) 

Qiiirkitit 

c Replace the input variables with a local variable so that 
c equations in the form of y=x+y will not lead to confusion, 
c 

zt1x=z1x 

zt2x=z2x 

c 

clogzh=0.5d0*(zt 1x-zt2x) 
dxh=ztmp(1 ) 

ifCdxh .It. zero) then 
zx=zt2x 
dxh=-dxh 
else 

zx=zt1x 
end if 
^«*****«* 

c machine accuracy = 2**(-53) 

c 2**(27)=e**e2p27 

c 

if (dxh .ge. e2p27) then 
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50 
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58 

59 
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70 
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73 
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76 

77 

78 

79 

80 



return 

else 

zx=0.5d0*(zt1x+zt2x) 
dsum=cdexp( c I ogzh ) 
dsum=1 .dO/dsum+dsum 
if (cdabs(dsum) .gt. ctol) then 
zx=cdlog(dsum)+zx 
else 

Cancellation is serious. ImCclogzh] is close to pi/2 or -pi/2. 
yi=dnint(ztmp(2)/twopi )*2.d0 
ztmp(2)=ztmp(2)-pi*yi 
dyi=dpi*yi 

if (ztmp(2) .It. zero) then 
clogzh=-clogzh 
dyi=-dyi 
end if 

ztmp(2)=(ztmp(2)-hpi )-hdpi-dyi 
dsum=2.d0*clogzh*(cone+c16*clogzh*clogzh) 
if (dsum .eq. czero) then 

Note that a complete cancellation of two nonzero numbers of 
order one is considered to be as accurate as what is allowed 
by the machine and the algorithm. 
zx=cerrx+chpi+zx 
else 

dsum=cd log (dsum) 

if (stmp(l) .It. e2m54) stmpCI )=e2m54 
zx=dsum+chpi+zx 
end if 
end if 
return 
end if 
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APPENDIX B: SUBROUTINE FZEROX 



This Appendix includes the listing of the subroutine FZEROX which identifies 
the meshes which may contain modes within a contour rectangle. The Shellman- 
Morffit mode locating algorithm has been completely replaced. 
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5 

6 

7 
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9 

10 

11 

12 

13 

14 

15 

16 
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18 

19 
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21 
22 

23 

24 

25 

26 

27 

28 

29 
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31 

1 

2 

3 
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32 

33 

34 

35 

36 

37 
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subroutine fzeroxCt left, t right, tbot, t top, tmshO, zeros, ni ,nf) 

kitit* 

fzerox is a routine for finding the zeroes of a complex function, f, 
which lie within a specified rectangular region of the 
complex q11 plane, assuming that the function has only 
simple zeroes over this rectangle. 

parameters specifying the search rectangle: 

tleft - value of the real part of q11 at the left edge, 

tright' value of the real part of q11 at the right edge, 

tbot - value of the imaginary part of q11 at the bottom edge, 
(this is set to 0. ) 

ttop • value of the inf^aginary part of q11 at the top edge, 

tmesh - set equal to about half the average spacing between 

zeroes within the rectangle. A smaller value may be used 
as a safety measure, but too small a value will result 
in excessively long run time, 
zeros - output list of (complex) values of q11 at which 
zeroes are found. 

nf-ni “ the number of zeroes found 

subroutines calledd-* 
f i ndfx 
roots 
nomshx 

itirkifk 



implicit double precision (a-h,o-z) 

comp I ex* 16 f 10, fOI , f 1 1 , fxnew, fxold, fxOO, fx10,fx01,fx11, 

+ czero, one, ci , sol , zeros 

parameter(czero=(0.d0,0.d0),one=(1 .d0,0.d0),ci=(0.d0, 1 .dO)) 
$ include: 'mlaparm. inc' 

***** Begin listing of: mlaparm. inc 
c 



include file to define the 

maximum # of layers (mxlayr) 
maximum # of modes (mxmode) 

parameter (mxlayr=35 ) 
parameter (mxmode=127) 

End listing of: mlaparm. inc 

dimensi on kedgel ( 100) ,kedge2( 100), kedge3( 100),kedge4( 100) , 

loci 2r( mxmode) , loc12i (mxmode), loc23r(mxmode), loc23i (mxmode), 
+ I oc34r( mxmode), loc34i (mxmode), I oc41r (mxmode), loc41 i (mxmode), 
sol(3),theta(2),zeros(2*mxmode+1 ) 



common /tmccom/ tmesh 
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59 
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61 
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64 

65 

66 

67 

68 

69 

70 

71 

72 

73 

74 

75 

76 

77 

78 

79 

80 

81 

82 

83 

84 

85 



^****«r 

c maxnsq - maximum number of mesh squares allowed on any one 
c phase line 

c maxnt - maximum number of times fzerox will reduce tmesh 
c 

maxnsq=3*max0( int ((ttop-tbot )/tmsh0), int ((tright- t lef t )/tmsh0)) 
maxnt=2 

tmesh = tmshO 
ntime = 0 
go to 7 
c 

5 tmesh=tmesh/2.0d0 

ntime = ntime+1 
if(ntime .gt. maxnt) go to 97 
c 

7 continue 

c calculate coordinates of rectangle edges in tmesh units 

c 

jit = idnintC t left/tmesh-0.5d0) 
jrt = idnintC tright/tmesh+0.5d0) 
jtop = idnintC ttop/tmesh+1 .5d0) 
jbot = 0 
c 

c initialize parameters for starting search at upper left 

c corner of search rectangle 

c 

ki = jtop 
kr = jit 
kedge = 1 

cal I f indfxCkr, ki , fxnew,xnew,ynew) 

nre1=0 

nre2=0 

nre3=0 

nre4=0 

knot 12=0 

knot23=0 

knot34=0 

knot41=0 

nf=ni 

n i 1 =n 1 + 1 

go to 15 

10 continue 

ifCnrzl .It. 2) go to 15 
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86 


c 


wri te(16,2000) nrzl 


87 




go to 5 


88 


15 


nrzl=0 


89 




nrsqu = 0 


90 


20 


fxold=fxnew 


91 




xold=xnew 


92 




yold=ynew 


93 




go to (21 ,26,31 ,36),kedge 


94 




95 


c 


search along left edge of rectangle for changes in the 


96 


c 


sign of imag(f) 


97 


c 




98 


21 


continue 


99 




if (ki ,eq. jbot) then 


100 




kedge=2 


101 




go to 26 


102 




end if 


103 




II 


104 




call f indfx(kr,ki , fxnew,xnew,ynew) 


105 




if (yold*ynew ,gt. O.dO) go to 20 


106 




if (nrel .eq.O) go to 23 


107 


c 




108 


c 


check if crossing point has been previously found 


109 


c 




no 




do 22 k=1 ,nre1 


111 




if (ki .eq.kedgel(k)) go to 20 


112 


22 


cont i nue 


113 


c 




114 


c 


follow phase line through rectangular region 


115 


c 




116 


23 


fx01=fxold 


117 




fxOI r=xold 


118 




fxOI i=yold 


119 




fx00=fxnew 


120 




fx00r=xnew 


121 




fx00i=ynew 


122 




li = ki 


123 




Ir = jit 


124 




go to 43 


125 




126 


c 


search along bottom edge of rectangle for changes in the 


127 


c 


sign of imag(f) 


128 


c 




129 


26 


continue 


130 




if (kr.eq. jrt) then 


131 




kedge=3 


132 




go to 31 
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133 

134 

135 

136 

137 

138 

139 

140 

141 

142 

143 

144 

145 

146 

147 

148 

149 

150 

151 

152 

153 

154 

155 

156 

157 

158 

159 

160 

161 

162 

163 

164 

165 

166 

167 

168 

169 

170 

171 

172 

173 

174 

175 

176 

177 

178 

179 



27 
c 
c 
c 

28 



CW1 

c 

c 

31 



32 
c 
c 
c 

33 



end if 
kr = kr+1 

cal I f indfx(kr,ki , fxnew,xnew,ynew) 
if (yold*ynew .gt. O.dO) go to 20 
if (nre2.eq.O) go to 28 

check if crossing point has been previously found 

do 27 k=1,nre2 
if (kr.eq.kedge2(k)) go to 20 
continue 

follow phase line through rectangular region 

fx00=fxold 
fx00r=xold 
fx00i=yold 
fx10=fxnew 
fx10r=xnew 
fx10i=ynew 
li = jbot 
Ir = kr-1 
go to 48 

**★ 

search along right edge of rectangle for sign changes in imag(f). 
continue 

if (ki ,eq. jtop) then 
kedge=4 
go to 36 
end if 
ki = ki+1 

cal I f indfx(kr , ki , f xnew,xnew,ynew) 
if (yold*ynew .gt. O.dO) go to 20 
i f (nre3.eq.O) go to 33 

check if crossing point has been previously found 

do 32 k=1,nre3 

if (ki .eq.kedge3(k)) go to 20 

continue 

follow phase line through rectangular region 

fx10=fxold 
fx10r=xold 
fx10i=yold 
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180 

181 

182 

183 

184 

185 

186 

187 

188 

189 

190 

191 

192 

193 

194 

195 

196 

197 

198 

199 

200 

201 

202 

203 

204 

205 

206 

207 

208 

209 

210 

211 

212 

213 

214 

215 

216 

217 

218 

219 

220 

221 

222 

223 

224 

225 

226 



fxl 1=fxnew 
fxl 1 r=xnew 
fxl 1 i=ynew 
li =: kf-1 
Ir = jrt“1 
go to 53 

c search along top edge of rectangle for sign changes in imag(f). 
c 

36 continue 
ifCkr.eq.jlt) go to 80 
kr = kr-1 

call f indfxCkr, ki , fxnew,xnew,ynew) 
if (yold*ynew .gt. O.dO) go to 20 
if (nre4.eq.O) go to 38 
c 

c check if crossing point has been previously found 

c 

do 37 k=l,nre4 
i f (kr.eq.kedge4(k)) go to 20 

37 continue 
c 

c follow phase line through rectangular region 

c 

38 fxl1=fxold 
fxllr=xold 
fxll i=yold 
fx01=fxnew 
fx01r=xnew 
fxOl i=ynew 
li = jtop-1 
Ir = kr 

go to 58 

c enter mesh square from left side or exit rectangle at right edge. 

41 lr=lr+l 

if (Ir .le. jrt-1) go to 42 
nre3=nre3+1 
kedge3(nre3)=l i+1 
go to 10 

42 fx01=fxll 
fx01r=fxl1r 
fx01i=fxl1i 
fx00=fxl0 
fx00r=fx10r 
fx00i=fx10i 
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227 

228 

229 

230 

231 

232 

233 

234 

235 

236 

237 

238 

239 

240 

241 

242 

243 

244 

245 

246 

247 

248 

249 

250 

251 

252 

253 

254 

255 

256 

257 

258 

259 

260 

261 

262 

263 

264 

265 

266 

267 

268 

269 

270 

271 

272 

273 



43 continue 

call f indfxC lr+1 , I i+1 , fxl 1 , fxl 1 r, fxl 1 i ) 
call f indfxC lr+1 , I i , f xIO, f xIOr, fxIOi ) 

^******4r 

c Determine the edge of exit of im(f)=0 from current mesh, 
edgei t = fx01 i’^fxl 1 i 
edgeib=fx00i’^fx10i 
if (edgeib ,gt. O.dO) then 
c lm(f)=0 goes through the 01 to 10 line, 

if (edgeit .gt. O.dO) then 

c lm(f)=0 goes through the 10 to 11 edge (edge 1). 

lout=1 
else 

c lm(f)=0 goes through the 01 to 11 edge (edge 2) 

lout=2 
end if 
else 

c lm(f)=0 goes through the 00 to 10 edge (edge 4) 

lout=4 

if (edgeit .It. O.dO) then 

c lm(f)=0 also runs through 01 to 11 and 10 to 11 edges, 

c Store crossing location and in/out information. 

knot34=knot34+1 
c loc34r(knot34)=lr 

c loc34i (knot34)=l i 

end if 
end if 
^******* 

go to 60 

Q-k’k'tt’k'k 

c enter mesh square from bottom side or exit rectangle at top edge. 

46 li=li+1 

if (li .le. jtop-1) go to 47 
nre4=nre4+1 
kedge4(nre4)=lr 
go to 10 

47 fx00=fx01 
fx00r=fx01r 
fx00i=fx01i 
fx10=fx11 
fx10r=fx11r 
fx10i=fx11i 

48 continue 

call f indfx( Ir , I i+1 , fxOI ,fx01 r ,fx01 i ) 
call f indfx( lr+1 , I i+1 , f x1 1 , fxl 1 r , fxl 1 i ) 

^******* 

c Determine the edge of exit of im(f)=0 from current mesh. 
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274 

275 

276 

277 

278 

279 

280 

281 

282 

283 

284 

285 

286 

287 

288 

289 

290 

291 

292 

293 

294 

295 

296 

297 

298 

299 

300 

301 

302 

303 

304 

305 

306 

307 

308 

309 

310 

311 

312 

313 

314 

315 

316 

317 

318 

319 

320 



edge! I=fx00i*fx01 i 
edge! r=fx10i*fx11 1 
if (edgeir .gt. O.dO) then 
c lm(f)=0 goes through the 00 to 11 line, 

if (edgeil .gt. O.dO) then 

c lm(f)=0 goes through the 01 to 11 edge (edge 2) 

lout=2 
else 

c lm(f)=0 goes through the 00 to 01 edge (edge 3). 

I out =3 
end if 
else 

c lm(f)=0 goes through the 10 to 11 edge (edge 1) 

lout=1 

if (edgeil .It. O.dO) then 

c lm(f)=0 also runs through 00 to 01 and 01 to 11 edges, 

c Store crossing location and in/out information. 

knot41=knot41+1 
c loc41r(knot41 )=lr 

c loc41 i (knot41 )=l i 

end i f 
end if 

Q******* 

go to 60 

c enter mesh square from right side or exit rectangle at left edge. 

51 lr=lr-1 

if (Ir .ge. jit) go to 52 
nre1=nre1+1 
kedgeKnrel )=l i 
go to 10 

52 fx11=fx01 
fx11r=fx01r 
fx11i=fx01 i 
fx10=fx00 
fx10r=fx00r 
fx10i=fx00i 

53 continue 

cal I f indfx( Ir, I i + 1 , fxOI , fxOI r , fxOI i ) 
call f indfxdr, I i , fxOO, fxOOr , fxOOi ) 

c Determine the edge of exit of im(f)=0 from current mesh. 
edgeit=fx01i*fx11i 
edgeib=fx00i*fx10i 
if (edgeit .gt. O.dO) then 
c lm(f)=0 goes through the 01 to 10 line. 
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321 

322 

323 

324 

325 

326 

327 

328 

329 

330 

331 

332 

333 

334 

335 

336 

337 

338 

339 

340 

341 

342 

343 

344 

345 

346 

347 

348 

349 

350 

351 

352 

353 

354 

355 

356 

357 

358 

359 

360 

361 

362 

363 

364 

365 

366 

367 



if (edgeib .gt. O.dO) then 

c lm(f)=0 goes through the 00 to 01 edge (edge 3). 

lout=3 
else 

c lm(f)=0 goes through the 00 to 10 edge (edge 4) 

lout=4 
end if 
else 

c lm(f)=0 goes through the 01 to 11 edge (edge 2) 

lout=2 

if (edgeib .It. O.dO) then 

c lm(f)=0 also runs through 00 to 10 and 00 to 01 edges, 

c Store crossing location and in/out information. 

knot12=knot12+1 
c loc12r(knot12)=lr 

c loc12i (knot12)=l i 

end if 
end i f 

go to 60 

c enter mesh square from top side or exit rectangle at bottom edge. 

56 li=li-1 

if ( I i .ge. jbot) go to 57 
nre2=nre2+1 
kedge2(nre2)=lr+1 
go to 10 

57 fx01=fx00 
fxOI r=fx00r 
fxOI i=fx00i 
fx11=fx10 
fx11r=fx10r 
fx11i=fx10i 

58 continue 

call f indfx( I r, I i , f xOO, fxOOr, fxOOi ) 
call findfx(lr+1,li,fx10,fx10r,fx10i) 

^4r****** 

c Determine the edge of exit of im(f)=0 from current mesh, 
edgei I=fx00i*fx01 i 
edgei r=f x10i*fx1 1 i 
if (edgei I .gt. O.dO) then 
c lm(f)=0 goes through the 00 to 11 line, 

if (edgeir .gt. O.dO) then 

c lm(f)=0 goes through the 00 to 10 edge (edge 4) 

lout=4 
else 

c lm(f)=0 goes through the 10 to 11 edge (edge 1). 
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368 

369 

370 

371 

372 

373 

374 

375 

376 

377 

378 

379 

380 

381 

382 

383 

384 

385 

386 

387 

388 

389 

390 

391 

392 

393 

394 

395 

396 

397 

398 

399 

400 

401 

402 

403 

404 

405 

406 

407 

408 

409 

410 

411 

412 

413 

414 



lout=1 
end if 
else 

c lm(f)=0 goes through the 00 to 01 edge (edge 3) 

lout=3 

if (edgeir ,lt. O.dO) then 

c lm(f)=0 also runs through 00 to 10 and 10 to 11 edges, 

c Store crossing location and in/ out information. 

knot23=knot23+1 
c loc23r(knot23)=lr 

c loc23i (knot23)=l i 

end if 
end if 
c 

60 continue 

nrsqu=nrsqu+1 

ifCnrsqu .gt. maxnsq) go to 95 

c Test for there being at least one re(f)=0 line entering and 
c leaving the mesh square, 

c 

if ((fx00r’*'fx10r .gt. O.dO) .and. (fx01r*fx11r .gt. O.dO) 

+ .and. (fxOOr^fxOIr .gt. O.dO)) go to (41,46,51,56) lout 
c 

c Computate the values of the modal function at the corners of a 
c a mesh square to determine its Taylor series to the 3rd order 

c for estimating its root locations, 

c 

c f00=one 

f 10=cdexp(fx10-fx00)-one 
f01=cdexp(fx01 -fx00)-one 
f 11=cdexp(fx1 1-fx00)-one 
c 

c write (16,3001) ni ,nf , Ir, li ,knot12,knot23,knot34,knot41 

c 3001 format (/' ni, nf, Ir, li and knot12, 23, 34 and 43 before ROOTS 
c + 2i6,2x,2i6,2x,4i6) 

c 

estimate locations of zeroes by radicals ********-*^*******-*^ 
c 

call roots(f 10,f01 ,f 11 ,sol ,nrsol ) 
c 

do 63 n=1,nrsol 

ureal = dreal (sol (n) ) 
uimag = dimag(sol (n) ) 

if (ureal .It. O.dO .or. ureal .gt. I.OdO) go to 63 
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A15 

A16 

417 

418 

419 

420 

421 

422 

423 

424 

425 

426 

427 

428 

429 

430 

431 

432 

433 

434 

435 

436 

437 

438 

439 

440 

441 

442 

443 

444 

445 

446 

447 

448 

449 

450 

451 

452 

453 

454 

455 

456 

457 

458 



if (uimag .It. O.dO .or. uimag .gt. l.OdO) go to 63 

62 thetad )=( I r+ureal )*tmesh 
theta(2)=( I i+uimag)*tmesh 
nf = nf+1 

zerosCnf )=dcmplx(theta(l ), theta(2)) 
nrzl=nrzl+l 

63 cont i nue 

^★★★★★★★★♦★♦★♦★★♦♦★♦★★* ********* ******* ********************** ★★*★*** 
c write (16,3002) ni,nf,nrsol 

c 3002 formate/' out of ROOTS at 63, ni, nf and # of roots ',3i4) 

Q* ****************************************************************** 

c continue following the phase line 

go to (41,46,51,56) lout 

Q****** 

cc 

80 cont i nue 

c 

return 

^***** 

95 continue 

wri te(16,9500) 

wri te(16,4001 ) I r , I i ,ni ,nf , tmesh 
write(* ,9500) 

4001 formate 'go to 5 from 95 at Ir, li = ' , i6, ' , ' , i6, ' ni, nf =',i6, 
+',',i6,', mesh size =',dl4.6) 

go to 5 

Q***** 

97 continue 

wri te(16,9700) 

write(16,4002)lr, I i ,ni ,nf, tmesh 
write(* ,9700) 

4002 formate 'go to 5 from 97 at Ir, li =' , i6, ' , ' , i6, * ni, nf =',i6, 
+',',i6, ', mesh size =' ,dl4. 6,/ ' zeroes found are kept.') 

c nf=ni 

c 

return 



c**** format statements 

9500 format e/5x, ' too many squares on same phase line -- ', 

$ 'reduce tmesh and start over') 

9700 formate/5x, ' tmesh has been reduced but problems remain in', 
$ ' executing fzerox') 

c 

end 
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APPENDIX C: SUBROUTINE ROOTS 



This Appendix contains the listing of the subroutine ROOTS. This subroutine 
replaces the portion of the subroutine FZEROX where the coefficients of a quadratic 
equation are determined, and the subroutine QUAD for locating the zeroes of a 
quadratic polynomial. In the revised subroutine FZEROX, the roots of a cubic 
polynomial has to be found. This subroutine determines these zeroes by radicals. 
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1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 

21 

22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 



subroutine roots f 1 , f 2, f3, sol , nrsol ) 

Qit it it it it it it it * -k ii It It if if -k -kirk it -k it it -kit it it it it -k It it h It hh it -k it ii-kii it it it it h -it -irirk it it it ■kit ■it it it -tr nitwit ■kit* it it * 

c This subroutine finds the roots of a third order polynomial by 
c radicals when the values of this polynomial at z=0, z=1, z=i and 
c z=Ui are given as f0=1, f1 + f0, f2+f0 and f3+f0 respectively, 
c Note that this algorithm takes cubic roots of two complex numbers 
c (hence the name 'solution by radicals') and use their linear 
c combinations as the roots of a third order polynomial. 

implicit real*8 (a-h, o-z) 

complex* 16 f 1 , f2, f3,zero,one,ci ,ep14,em14,ep23,em23, 

f a, fb, fc, fd, f a1 ,fa2,fa3, fa1s,p,q,delt,z,zm,u,v,sol 
parameter (xbi t52=52.d0*0.69314718055994531d0, thrd=1 .d0/3.d0, 

+ bit50=1.d0/33554432.d0/33554432.d0,bit51=bit50/2.d0, 

+ bit52=bit51/2.d0,tol=0.001d0, 

zero=(0.d0,0.d0),one=(1 .d0,0.d0),ci = (0.d0,1 .dO), 

+ ep14=(0.5d0,0.5d0),em14=(0.5d0,-0.5d0), 

+ ep23= ( - 0 . 5d0 , 0 . 86602540378443864675d0 ) , 

em23=( - 0 . 5d0 , - 0 . 86602540378443864675d0 ) ) 
dimension sol(*) 
f a=one 

fb=(f2-ci*f1+em14*f3) 
f c=((ep14+one)*f 1 - (em14+one)*f 2+ci*f3) 



fd= 


:(emU*(f2-f1)-epU*f3) 




if 


(cdabs(f b) 


.le. bit50) 


f b=zero 


if 


(cdabs(fc) 


.le. bit51) 


f c=zero 


if 


(cdabs(fd) 


.le. bit52) 


fd=zero 


if 


(fd ,ne. zero) then 






fa1=(- thrd) 


*fc/fd 





fa2=fb/fd 



fa3=fa/fd 
fa1s=fa1*fa1 
p=thrd*f a2-fa1s 

q=0.5d0*(fa3+fa1*fa2)-fa1*f als 
if (p .eq. zero) then 
if (q. eq. zero) then 
nrsol=1 
sol ( 1 )=fa1 
return 
else 

nrsol=3 

u=((-2.d0)*q)**thrd 
sold )=u+fa1 
sol (2)=ep23*u+fa1 
sol (3)=em23*u+f a1 
return 
end if 
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48 

49 

50 

51 

52 

53 

54 

55 

56 

57 

58 

59 

60 

61 

62 

63 

64 

65 

66 

67 

68 

69 

70 

71 

72 

73 

74 

75 

76 

77 

78 

79 

80 

81 

82 

83 

84 

85 

86 

87 

88 

89 

90 

91 

92 

93 

94 



else 



if (q. eq. zero) then 
nrsol=3 
sold ) = fa1 
u=cdsqrt( (-3 .d0)*p) 
sol (2)=fa1+u 
sol(3)=fa1 -u 
return 
else 
v=p/q 
Z=p* v*v 
absz=cdabs(z) 
if (absz .It. tol) then 
zm=-z 

fn=dint(1 .dO-xbi t52/dlog(absz) ) 

lastn=idint(fn)- 1 

dnn=fn-0.5d0 

dnd=fn+1.0d0 

delt=one 

do 100 nt=1 , lastn 
dnn=dnn- 1 .dO 
dnd=dndd.d0 

de 1 1 = ( dnn/ dnd ) * de 1 1 * zfTH one 
continue 

del t=(0.5d0*del t/q)**thrd 
u=p*delt 
v=-1 .dO/delt 
else 

de 1 1 =cdsqr t ( one+ z ) - one 
u=(q*delt )**thrd 
v=-p/u 
end if 
nrsol=3 

sol ( 1 )=u+v+fa1 
sol (2)=ep23*u+em23*v+fa1 
sol (3)=em23*u+ep23*v+fa1 
return 
end if 
end if 

else if (fc .ne. zero) then 
if (fb .eq. zero) then 
if (fa .eq. zero) then 
nrsol=1 
sol ( 1 )=zero 
return 
else 

nrsol=2 
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95 

96 

97 

98 

99 

100 

101 

102 

103 

104 

105 

106 

107 

108 

109 

110 

111 

112 

113 

114 

115 

116 

117 

118 

119 

120 

121 

122 

123 

124 

125 

126 

127 

129 

130 

131 

132 

133 

134 

135 

136 

137 

138 



z=cdsqrt(-fa/fc) 
sol ( 1 )=z 
sol(2)=-z 
return 
end if 
else 

fa1=0.5d0*fb/fc 

fa2=fa/fc 

Z=fa2/fa1/fa1 

absz=cdabs(z) 

if (absz .It. tol) then 

fn=dint( 1 .dO-xbi t52/dlog(absz) ) 

lastn=idint(fn)- 1 

dnn=fn-0.5d0 

dnd=fn+1 .OdO 

del t=one 

do 200 nt=1 , lastn 
dnn=dnn- 1 .dO 
dnd=dnd-1.d0 

del t=(dnn/dnd)*del t*z+one 
continue 

delt=-0.5d0*delt/fa1 

nrsol=2 

sole 1 ) = fa2*delt 
sol(2)=1 .dO/del t 
return 
else 

del t=cdsqr t ( one- z) 
nrsol=2 

sol ( l)=-fa1 *( one-del t) 
sol(2)=-fa1*(one+del t ) 
return 
end if 
end if 

else if (fb .ne. zero) then 
nrsol=1 
sol(1)=-fa/fb 
return 
else 

nrsol=1 
sole 1 )=ep14 
return 
end if 
end 
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APPENDIX D: SUBROUTINE ABCOEF 



This Appendix contains the listing of the subroutine ABCOEF. The consistency 
self-checking procedure has been implemented to determine the correct method to 
evaluate the /!,• and B,- coefficients. 
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1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 

21 

22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 



subroutine abcoef (zero,m) 



c For each mode m, this suboutine calculates A-B coefficients in 

c all layers for combining two linearly independent solutions of 

c Stokes' equation to form the height gain function: 
c 

c height gain=exp(bcoefx(l,m))*(k1*exp(acoefx(l,m))+k2) 

c 

c where k1 and k2 are two independent solutions to Stokes' 

c equation. In the top layer (i.e. nzlayr) the height gain is: 

c 

c height gain=exp(bcoefx( I ,m))*h2 

c where h2 is a solution to the Stokes' equation associated 
c with outgoing energy flow. Here k1 and k2 are proportional 

c to the k1 and k2 used by Marcus and the h2 is proportional 

c to a modified Hankie function of order 1/3. 

c inputs... 

c zero-an eigenvalue in q11 space 

c outputs... 

c acoefx-two dimensional array of complex exponents 

c coefficients used to combine two linearly 

c independent solutions of stokes' equation 

c bcoefx-two dimensional array of complex exponents 

c coefficients used for normalizing the height gains 

c note: acoefx and bcoefx are passed by the 

c common block /pap2/ 

c subroutines called... 

c xcdai 

c xcadd 

c common block areas... 

c comi 

c com2 

c papi 

c pap2 



implicit real*8(a-h,o-z) 

complex*16 acoefx, bcoefx, cqi j ,h2xq1 ,dh2xq1 ,h2xq2,dh2xq2, klxql , 
$ dklxql ,k1xq2,dk1xq2, k2xq1 ,dk2xq1 , k2xq2,dk2xq2,h2dk1x, 

$ dh2k1x,h2dk2x,dh2k2x,numax,denax,numbx,denbx, intlx, int2x, 

$ hyx,dhyx,k1dhyx,dk1hyx,dk2hyx,k2dhyx, gamma, dgamdq, i , 
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48 $ koa123, rtSLmx, zero,q1 ,q2,sumx, surfno,dqi j ,dqi jdz, sqng, 

49 $ dnumbx,dhux,dhlx,e13x,cneg,cldqzl ,cldqzm,cigama,koawav, tthd, 

50 + tacoef ,dacoef 

51 

52 parameter(downi=1 .d-3,downr=1 .d-3/0.4342944819032518d0, 

53 + pi =3 . 1 4 1 592653589793238462643d0 , 

54 + i=(0.0d0,1.0d0),tthd=(2.d0/3.d0)*i, 

55 + cneg=(0.0d0,3.141592653589793238462643d0),e13x=cneg/3.d0) 

56 c***** 

57 c mxlayr=maximum number of layers allowed 

58 c mxmode=maxfmLm number of modes allowed 

59 

60 c 

61 c use include file for parameters of 

62 c use include file for parameters of 

63 c mxlayr max # layers 

64 c mxmode max # modes 

65 c 

66 Sinclude: 'mlaparm. inc' 

***** Begin listing of: mlaparm. inc 

1 c 

2 c include file to define the 

3 c maximum # of layers (mxlayr) 

4 c maximum # of modes (mxmode) 

5 c 

6 parameter (mxlayr=35 ) 

7 parameter (mxmode=127) 

***** End listing of: mlaparm. inc 

67 c 

68 c 

69 c***** 

70 c acoefx-two dimensional complex array used for combining two 

71 c independent solutions to stokes' equation 

72 c bcoefx'two dimensional complex array used for normalizing height 

73 c gain 

74 c cqij-two dimensional array containing coefficients for evaluating 

75 c qij in terms of q11 

76 c dqij -array containing coefficients for evaluating qij in terms of 

77 c q11 

78 c dqi jdz-array containing derivatives of qi(z) in the different 

79 c layers 

80 c zi-array containing input hesights for the modified refractivity 

81 

82 dimension acoefx (mxlayr, mxmode), 

83 $ bcoefx (mxlayr, mxmode), 

84 $ dqi j (mxlayr), cqi j(mxlayr, 2) , dqi Jdz(mxlayr) , zi (mxlayr+1 ) 

85 c***** 
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87 

88 
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106 
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110 

111 
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113 
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118 
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125 

126 

127 

128 
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130 

131 

132 



common /com1/f req, waveno, sqng 
common /com2/cqi j,dqi j,dqi Jdz,nzlayr 
common /papl/nrmode, koal 23 , surf no, z i 
common /pap2/acoefx,bcoefx 



c check for single layer 

c 

c set a complex variable koawav=- i*koa123/(waveno*waveno) to 
c avoid repeating computations 

koawav= - i * koa 1 23/ ( waveno* wa veno ) 

if(nzlayr .eq. 1)then 
q1=cqi j(1 , 1)+zero*dqi j(1 ) 
call surf (q1 , gamma, dgamdq) 

cal I xcdai(-q1 ,k2xq1 ,dk2xq1 ,k1xq1 ,dk1xq1 ,h2xq1 ,dh2xq1 ) 
dh2xq1 =dh2xq1 +e1 3x 

i nt 1 x=cdl og( koawav*dgamdq- q1 /dqi jdz ( 1 ) ) +2 . 0d0*h2xq1 
int2x=2.0d0*dh2xq1 -cdlogC -dqi jdz( 1 ) ) 
call xcadd(sumx, int lx, int2x) 
rtsumx=0.5d0*sumx 
bcoefxC 1 ,m)=-rtsumx 
return 
end if 

cldqzl=cdlog(-dqi jdz(1)) 

c if I equals one then initialize cumulants and caculate a's and 
c b*s in bottom layer using ground boundary conditions. 

q1=cqi j(1 ,1)+zero*dqi j(1) 

cal I xcdai ( -q1 , k2xq1 ,dk2xq1 ,k1xq1 ,dk1xq1 ,h2xq1 ,dh2xq1 ) 

dk2xq1=dk2xq1+cneg 

dk1xq1=dk1xq1 -e13x 

call surf (q1 , gamma, dgamdq) 

cigama=cdlog( i*gamma) 

cal I xcadd(numax,cldqzl-cneg+dk2xq1 ,cigama+cneg+k2xq1 ) 

cal I xcadd(denax,cigama+k1xq1 ,cldqzl+dk1xq1 ) 

acoefxd ,m)=numax“denax 

cal I xcadd(denbx,k2xq1 ,acoefx(1 ,m)+k1xq1 ) 

bcoefxC 1 ,m)=-denbx 

c calculate contributions to normalizing integrals, 
cal I xcaddChyx, k2xq1 ,acoefx( 1 ,m)+k1xq1 ) 
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hyx=bcoef x( 1 ,m)+hyx 

call xcacid(dhyx,dk2xq1 ,acoefx(1 ,m)+dk1xq1 ) 
dhyx=bcoefx(1 ,m)+dhyx 

int1x=cdlog(koawav*dgamdq-q1/dq1 jdz(1) )+2.0d0*hyx 

int2x=2.0d0*dhyx*cldqzl 

call xcadd(sunx, intlx, int2x) 

do 9010 l=2,nzlayr-1 
lm1=l-1 

cldqzl=cdlog(-dqi jdz( D) 
cldqzm=cdlog(dqi jdz( Iml )) 
q1=cqi j(l,1)+zero*dqi j(l) 

call xcdaf <-q1 ,k2xq1 ,dk2xq1 ,k1xq1 ,dk1xq1 ,h2xq1 ,dh2xq1 ) 

dk2xq1=dk2xq1+cneg 

dk1xq1=dk1xq1-e13x 

q2=cqi J ( Iml , 2)+zero*dqi j ( Iml ) 

cal I xcdai C-q2,k2xq2,dk2xq2,k1xq2,dk1xq2,h2xq2,dh2xq2) 

dk2xq2=dk2xq2+cneg 

dk1xq2=dk1xq2*e13x 

cal I xcaddChyx, k2xq2,acoefx( Iml ,m)+k1xq2) 

call xcadd(dhyx,dk2xq2,acoefx( Iml ,m)+dk1xq2) 

k1dhyx=k1xq1+dhyx 

dk1hyx=dk1xq1+hyx 

dk2hyx=dk2xq1 +hyx 

k2dhyx=k2xq1 +dhyx 

cal I xcadd(denax,cldqzm+k1dhyx,cldqzl+dk1hyx) 

ca 1 1 xcaddC numax , c I dqz I - cneg+dk2hyx , c ldqzm+cneg+k2dhyx ) 

acoefxC I ,m)=numax-denax 

call xcadd(denbx,k2xq1 ,acoefx(l,m)+k1xq1 ) 

numbx=bcoefx( Iml ,m)+hyx 

dnumbx=bcoefx( Iml ,m)+dhyx 

bcoefxC I ,m)=numbx-denbx 

calculate contribution to normalizing integrals. 

int1x=cdlog(-q1/dqi jdz( I )+q2/dqi jdz( Iml ))+2.0d0*numbx 
call xcadd(sumx,sumx, intlx) 
call xcadd(dhux,dk2xq1 ,acoefx( I ,m)+dk1xq1 ) 
dhux=bcoef x( I , m)+dhux 
i nt 1 x=2 . 0d0*dnumbx - c I dqzm 
i nt 2x=2 . 0d0*dhux- c I dqz I 
call xcadd(sumx,sumx, intlx) 
call xcadd(sLinx,sumx, int2x) 
cont i nue 



if I equals nzlayer, calculate a's and b's using outgoing 
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c wave in top layer. 

nzm1=nzlayr- 1 

q1=cqi j (nzlayr, 1 )+zero*dqi j (nzlayr) 

call xcdai (-q1 ,k2xq1 ,dk2xq1 ,k1xq1 ,dk1xq1 ,h2xq1 ,dh2xq1 ) 

dh2xq1=dh2xq1+e13x 

q2=cqi j(nzm1 ,2)*»’Zero*dqi jCnzmI ) 

call Xcdai(-q2,k2xq2,dk2xq2,k1xq2,dk1xq2,h2xq2,dh2xq2) 

dk2xq2=dk2xq2+cneg 

dk1xq2=dk1xq2-e13x 

cal I xcadd(hyx,k2xq2,acoefx(nzm1 ,m)+k1xq2) 
numbx=bcoefx(nzlayr-1 ,m)+hyx 
bcoefx(nzlayr,m)=numbx-h2xq1 

c calculate contribution to cumulants. 

int1x=cdlog( -ql/dqi jdz(nzlayr)+q2/dqi jdzCnzmI ))+ 

$ 2.0d0*numbx 

call xcadd(sumx,sumx, intlx) 

cal I xcadd(dhyx,dk2xq2,acoefx(nzm1 ,m)+dk1xq2) 

dnumbx=bcoefx(nzm1 ,m)+dhyx 

int 1x=2.0d0*dnumbx-cdlog(dqi jdz(nzm1 ) ) 

call xcadd(sumx,sumx, intlx) 

dhux=bcoefx(nzlayr,m)*»'dh2xq1 

int2x=2.0d0*dhux-cdlog(-dqi jdz(nzlayr) ) 

ca 1 1 xcaddC sumx , sumx , i nt2x ) 

c renormalize b's so that height gain integral equals unity. 

r t s umx = . 5 d 0 s umx 
do 9000 1 1=1 , nzlayr 

bcoefxC 1 1 ,m)=bcoefx( 1 1 ,m)-rtsumx 
9000 continue 



l=nzlayr 

lm1=l-1 

cldqzm=cdlog(dqi jdzClml ) ) 
cldqzl=cdlog( -dqi Jdz( I )) 

c calculate q and associated quantities at bottom of layer I 
q1=cqi j(l ,1)+zero*dqi j(l) 

cal I xcdai (-q1 ,k2xq1 ,dk2xq1 , k1xq1 ,dk1xq1 ,h2xq1 ,dh2xq1 ) 
dh2xq1=dh2xq1+e13x 



61 



227 

228 

229 

230 

231 

232 

233 

234 

235 

236 

237 

238 

239 

240 

241 

242 

243 

244 

245 

246 

247 

248 

249 

250 

251 

252 

253 

254 

255 

256 

257 

258 

259 

260 
261 
262 

263 

264 

265 

266 

267 

268 

269 

270 

271 

272 

273 



q2=cqi j ( Iml ,2)+zero*dqi j ( Iml ) 

call xcdai ( *q2, k2xq2,dk2xq2,klxq2,dklxq2,h2xq2,dh2xq2) 

dk2xq2=dk2xq2+cneg 

dklxq2=dklxq2-el3x 



Q* ** ** 

c Caculate acoefxC Iml ,m),bcoefx( Iml ,m) 

c and cumulants using outgoing wave in nzlayr 

Q**irk* 

dh2klx=dh2xql+klxq2 

h2dklx=h2xql+dklxq2 

h2dk2x=h2xql+dk2xq2 

dh2k2x=dh2xql+k2xq2 

call xcaddC denax , c I dqz I - cneg+dh2k 1 x , c I dqzfrHcneg+h2dk 1 x ) 
call xcadd(numax,cldqzrrH-h2dk2x,cldqzl+dh2k2x) 

c If in the nzlayr-1 layer the magnitudes of A coefficients from 
c integration up and down differ by less than 0.02 dB and their 
c phases differ by less than O.OOlpi, the A and B coefficients 
c obtained from integration up will be accepted. 

tacoef=numax-denax 
dacoef=tacoef -acoefxC Iml ,m) 
difr=dabs(dreal(dacoef )) 
if (difr .It. downr) then 
dif i=dimag(dacoef )/pi 
di f i=dabs(di f i -dnintCdi f i/2.dO)*2.dO) 
if (difi .It. downi) return 
end i f 

acoefxC Iml ,m)=tacoef 

call xcaddCdenbx,k2xq2, acoefxC Iml ,m)+klxq2) 
bcoefxC Iml ,m)=h2xql -denbx 

c calculate contributions to cumulants 

sumx=cdlogC-ql/dqi JdzC I )+q2/dqi jdzC Iml ) )+2.0d0*h2xq1 

call xcaddCdh I x,dk2xq2, acoefxC Iml ,m)+dklxq2) 

dhlx=bcoefxC Iml ,m)+dhlx 

intlx=2.0d0*dh2xql-cldqzl 

call xcaddC intlx,sumx, intlx) 

int2x=2.0d0*dhlx-cldqzm 

call xcaddC sumx, intlx, int2x) 

do 9030 l=nzlayr-l,2,-l 



62 



274 

275 

276 

277 

278 

279 

280 

281 

282 

283 

284 

285 

286 

287 

288 

289 

290 

291 

292 

293 

294 

295 

296 

297 

298 

299 

300 

301 

302 

303 

304 

305 

306 

307 

308 

309 

310 

311 

312 

313 

314 

315 

316 

317 

318 

319 

320 



lm1=l'1 

clciqzl=cdlog(-dqi jdz( D) 
cldqzm=cdlog(dqi jdz( Iml ) ) 

c calculate q and associated quantities at bottom of layer I 
q1=cqi j( I, 1 )+zero*dqi j( I ) 

cal I xcdai (-q1 ,k2xq1 ,dk2xq1 ,k1xq1 ,dk1xq1 ,h2xq1 ,dh2xq1 ) 

dk2xq1=dk2xq1+cneg 

dk1xq1=dk1xq1-e13x 

q2=cqi j ( Iml , 2)+zero*dqi j ( Iml ) 

cal I xcdai (-q2,k2xq2,dk2xq2,k1xq2,dk1xq2,h2xq2,dh2xq2) 

dk2xq2=dk2xq2+cneg 

dk1xq2=dk1xq2*e13x 

dh2xq2=dh2xq2+e13x 



c Calculate acoefxC Iml ,m) ,bcoefx( Iml ,m) and cumulants 

c using continuity relations in terms of the linearly 

c independent functions k1 and k2 

call xcadd(hyx,k2xq1 ,acoefx( I ,m)+k1xq1) 

call xcadd(dhyx,dk2xq1 ,acoefx( I ,m)+dk1xql ) 

k1dhyx=k1xq2+dhyx 

dk1hyx=dk1xq2+hyx 

dk2hyx=dk2xq2+hyx 

k2dhyx=k2xq2+dhyx 

call xcadd(denax,cldqzl -cneg+k1dhyx,cldqzfTH-cneg+dk1hyx) 

cal I xcaddCnumax, cldqznrH'dk2hyx,cldqzl + k2dhyx) 

acoef x( Iml ,m)=numax-denax 

call xcadd(denbx, k2xq2,acoefx( Iml ,m)+k1xq2) 

numbx=bcoef x( I ,m)+hyx 

dnumbx=bcoefx( I ,m)+dhyx 

bcoefx( Iml ,m)=numbx-denbx 

c calculate contributions to cumulants. 

int1x=cdlog(-q1/dqi Jdz( I )+q2/dqi jdz( Iml ) )+2.0d0*numbx 

call xcadd(SLmx,sumx, intlx) 

cal I xcadd(dhlx,dk2xq2,acoefx( Iml ,m)+dk1xq2) 

dhlx=bcoefx( Iml ,m)+dhlx 

i n 1 1 x=2 . OdO^dnumbx - c I dqz I 

int2x=2.0d0*dhlx-cldqzm 

call xcadd(sunx,sLinx, intlx) 

call xcadd(sunx,sumx, int2x) 
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9030 



continue 



c if I equal to one calculate ground 
c contribution to cumulants and renormalize bcoefx's 

1 = 1 

q1=cqi j( I , 1 )+zero*dqi j( 1) 

call xcdai (-q1 ,k2xq1 ,dk2xq1 , klxql ,dk1xq1 ,h2xq1 ,dh2xq1 ) 

dk2xq1=dk2xq1+cneg 

dk1xq1=dk1xq1-e13x 

call xcaddChyx, k2xq1 , acoefxC I ,m)+k1xq1 ) 
cal I xcadd(dhyx,dk2xq1 ,acoefx( I ,m)+dk1xq1 ) 
call surf (q1 , gamma, dgamdq) 
numbx=bcoef x( I ,m)+hyx 
dnumbx=bcoefx( I ,m)+dhyx 

int1x=cdlog(koawav*dgamdq-q1/dqi jdz( I ) )+2.0d0*numbx 
int2x=2.0d0*dnumbx-cdlog(-dqi Jdz(1)) 
call xcadd(sumx,sumx, intlx) 
call xcadd(sumx,sumx, int2x) 

c renormalize b's so that height gain integrals equal unity. 

r t s umx = . 5 d 0 * s umx 

do 9020 1 1 = 1 ,nzlayr-1 

bcoefxC 1 1 ,m)=bcoefx( 1 1 ,m)-rtsumx 
9020 continue 

bcoefx(nzlayr,m)=-rtsumx 



return 

end 
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